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ABSTRACT 

In  this  thesis  we  study  several  computer  implementa- 
tions of  the  thinning  algorithm,  a new  method  for  generating 
non- homogeneous  Poisson  processes.  The  method,  developed 
by  Professor  P.A.W.  Lewis,  Naval  Postgraduate  School, 
Monterey,  California,  and  G.S.  Shedler,  IBM  Research 
Laboratory,  San  Jose,  California,  is  valid  for  Poisson 
processes  with  any  given  intensity  function.  The  basic 
thinning  algorithm  is  modified  to  exploit  several  refine- 
ments which  reduce  computer  execution  time  by  approximately 
one-third.  The  basic  and  modified  thinning  programs  are 
compared  with  a previous  algorithm  of  Lewis  and  Shedler, 
the  Poisson  decomposition  and  gap-statistics  algorithm, 

which  is  easily  implemented  for  Poisson  processes  with 

. . 2 
intensity  functions  of  the  form  exp (aQ+a1t+a2t  ).  The 

thinning  programs  are  competitive  in  both  execution  time 
and  computer  memory  requirements.  One  program  implementa- 
tion generates  the  events  in  a Poisson  process  one  at  a 
time;  another  program  implements  the  algorithmic  refinements 
which  improve  efficiency. 
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I . INTRODUCTION 


The  Poisson  process  is  a widely  known  and  studied 
stochastic  process.  It  is  frequently  used  to  describe 
random  arrivals  at  some  type  of  service  facility  such  as 
a service  station  fuel  pump  or  a bank  teller's  window. 

In  its  most  common  form,  the  "rate"  of  these  arrivals  is 
considered  to  be  constant  over  time.  This  is  the  homogeneous 
Poisson  process  which  has  the  familiar  property  that  times 
between  arrivals  (or  events)  are  exponentially  distributed 
with  mean  equal  to  the  inverse  of  the  rate. 

The  assumption  of  a constant  rate,  or  homogeneity,  is 
at  best  tenuous  when  applied  to  real  world  data.  For  exam- 
ple, the  rate  of  arrivals  at  a traffic  light  typically 
varies  from  very  high  during  rush  hours  to  very  low  in  the 
early  morning.  In  addition  to  this  cyclic  time-of-day 
effect,  arrival  rates  may  exhibit  longer  term  increases  or 
decreases.  Further,  these  effects  may  be  superimposed  upon 
shorter  term  effects  to  produce  a more  complex  rate  which 
varies  with  time.  These  processes  for  which  arrival  rates 
vary  with  time  may  often  be  represented  by  a non-homogeneous 
Poisson  process,  that  is.,  a Poisson  process  with  a time 
dependent  rate  of  arrival. 

The  generic  term  of  Poisson  process  includes  then  both 
homogeneous  and  non-homogeneous  Poisson  processes.  LEWIS 
and  SHEDLER  [Ref.  1]  define  the  Poisson  process  generally  in 
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terms  of  a monotone  non-decreasing  right-continous  function 

A(t)  which  is  bounded  in  any  finite  interval.  Then  the 

number  of  points,  N(t,,,t')/  in  any  finite  interval  has  a 

Poisson  distribution  with  parameter  u(t'',t')  = A(t')  - A(t''). 

Thus,  for  example  in  ( 0 , t ' ] , with  t'  >_  0,  P{N(t'',t')  = n} 

_yn 

= P{Nt,=n}  = uQ  e /n! , where  = y(0,t')  = A(t')  - A(0). 

The  right  derivative  A (t)  of  A(t)  will  be  assumed  to 
exist  and  is  called  the  rate  function  or  intensity  function 
of  the  process.  A(t)  is  called  the  integrated  rate  function 
and  has  the  interpretation  that  for  t ^ 0,  A(t)  - A(0)  = E [N^  ] . 
For  the  homogeneous  Poisson  process  A(t)  is  a constant, 
e.g.  X,  and  thus  the  integrated  rate  function  is  simply  the 
product  of  A and  t,  i.e.  the  expected  value  of  = N(0,t). 

While  simulation  of  homogeneous  Poisson  processes  is 
relatively  straightforward,  the  non-homogeneous  Poisson 
process  is  more  problematical.  Times  between  events  are 
not  exponential  in  the  general  case  and  simulation  has 
typically  been  tailored  to  specific  classes  of  intensity 
functions.  LEWIS  and  SHEDLER  [Ref.  1]  list  three  general 
methods  for  simulating  non-homogeneous  Poisson  processes 
and  one  method  for  a special  rate  function.  The  general 
methods  include  the  time  scale  transformation  method  and 
the  conditioning  and  order-statistics  method.  The  special 
method  is  the  gap-statistics  method,  a method  which  is 
particular  to  the  degree-one  exponential  polynomial  intensity 
fucntion,  i.e.  those  of  the  form  A (t)  = exp(bQ  + b^t)  . 
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Implementation  of  the  general  methods  on  a computer  may 
pose  special  problems . Often  the  inverse  of  the  integrated 
rate  function  is  not  explicit  and  must  be  computed  numerically. 
Other  problems  in  implementation  generally  result  in  lower 
efficiency,  as  measured  by  execution  time  or  computer  storage 
requirements  or  both. 

One  class  of  intensity  functions  which  is  of  general 

interest  is  the  degree-two  exponential  polynomial  family. 

That  is,  those  with  intensity  function  of  the  form 

2 

A(t)  = exp(aQ  + a^t  + a.^t  ).  This  family  of  functions  has 
the  property  of  being  positive  for  all  values  of  t,  a 
necessary  condition  for  an  intensity  function.  Additionally, 
by  varying  the  magnitude  and  sign  of  the  coefficients,  the 
exponential  polynomial  of  degree  two  can  be  made  to  be  mono- 
tone increasing  or  decreasing  over  time,  as  well  as  increasing 
and  then  decreasing,  or  vice  versa.  Use  of  this  type  of 
intensity  function  also  leads  to  statistical  procedures 
which  are  relatively  simple. 

LEWIS  and  SHEDLER  [Ref.  2]  proposed  a new  method  of 
generating  the  non-homogeneous  Poisson  process  with  degree- 
two  exponential  polynomial  intensity  function.  It  involves 
decomposition  of  the  degree-two  exponential  polynomial 
intensity  function,  A(t),  into  two  functions,  a degree-one 
exponential  polynomial  function,  AT  (t) , and  a difference 
function,  XD ( t ) = \ (t)  - AL(t).  This  procedure  allows  the 
points  in  the  degree-one  exponential  polynomial  event  stream 


to  be  generated  using  the  gap-statistics  method,  which  is 
highly  efficient  when  implemented  on  a computer.  The 
remaining  points  with  intensity  function  XD(t)  are  generated 
by  other  methods  and  then  merged  with  the  other  events. 

PATROW  [Ref.  3]  implemented  two  algorithms,  the  time 
scale  transformation  algorithm  and  the  Poisson  decomposi- 
tion and  gap-statistics  technique,  and  compared  them  for 
computational  speed  and  computer  memory  requirements.  His 
results  indicated  that  the  Poisson  decomposition  and  gap- 
statistics  technique  was  from  two  to  seven  times  faster  than 
the  time  scale  transformation  algorithm,  although  the  former 
required  about  thirty  percent  more  computer  memory. 

PATROW' s work  [Ref.  3]  is  also  an  excellent  self- 
contained  reference  on  Poisson  processes,  bringing  many 
references  together  under  one  cover. 

A recent  result  of  LEWIS  and  SHEDLER  [Ref.  1]  develops 
a new  method  for  generation  of  points  in  a non-homogeneous 
Poisson  process.  This  method,  called  "thinning",  is  similar 
to  the  general  conditioning-acceptance-rejection  method 
but  has  subtle  differences  which  are  computationally  signi- 
ficant. The  thinning  method  is  straightforward  in  both  an 
analytical  and  a computational  sense,  and  is  valid  for  any 
type  of  intensity  function.  The  thinning  theorem  is 
presented  in  Section  II. 

This  thesis  is,  in  a sense,  a sequel  to  PATROW 's  work 
[Ref.  3].  Its  purpose  is  to  implement  the  thinning  algorithm 
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in  computer  program  form  and  to  compare  it  to  the  Poisson 
decomposition  and  gap-statistics  algorithm  implemented  by 
PATROW  [Ref.  3].  The  latter  implementation  was  designed 
for  a specific  subset  of  intensity  functions,  degree-two 
exponential  polynomials.  Since  the  Poisson  decomposition 
and  gap-statistics  method  outperformed  a general  case 
algorithm  (time  scale  transformation)  by  a significant 
margin,  comparing  the  thinning  method  to  the  Poisson 
decomposition  and  gap-statistics  method  should  give  a 
reasonable  indication  of  the  thinning  algorithm's  performance 
in  generating  non-homogeneous  Poisson  processes  with  other 
than  degree-two  exponential  polynomial  intensity  functions. 

Section  III  lists  the  two  algorithms  considered,  as  well 
as  a special  application  of  the  thinning  process  which  will 
be  of  interest  to  those  involved  in  event-step  simulation. 

Section  IV  describes  the  methodology  used  in  comparing  the 
algorithms  while  Section  V deals  with  aspects  of  the  thinning 
procedure  which  may  be  exploited  to  enhance  its  overall 
effectiveness  in  a variety  of  situations.  Finally,  Section 
VI  presents  the  results  and  conclusions  of  the  comparisons 
of  the  algorithms.  Appendices  A and  B contain  secondary 
results  and  computer  program  listings  following  the 
appendices . 

i| 
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II.  THE  THINNING  THEOREM 


The  underlying  concept  of  the  thinning  method  involves 

★ 

the  use  of  a "bounding"  Poisson  process,  {Nt  : t >_  0}  , 

★ 

where  Nfc  is  the  number  of  points  in  the  bounding  process  in 

the  interval  (0,t].  This  process  may  be  either  homogeneous 

or non-homogeneous  Poisson,  but  should  be  one  which  is  easy 

to  simulate  on  a computer.  It  is  called  bounding  because 

★ 

its  intensity  function,  denoted  A (t) , bounds  the  intensity 
function  A(t),  of  the  nonhomogeneous  Poisson  process  which 

is  to  be  simulated  over  the  fixed  interval  ( 0 , t * ] . That 

* * 

is,  X (t)  A(t)  for  all  t in  (0,t']«  Points  at  T^ , 

★ 

i = 1,  ...,  Nt,,  are  generated  for  the  bounding  process 

over  the  interval  ( 0 , t * ] . These  points  are  then  deleted, 

or  "thinned",  with  independent  probabilities  equal  to 
★ * * 

1 - (X(T^)/A  (T^) ) . Thus  the  probability  that  a point  of 

★ 

the  bounding  process,  T^ , is  a point  of  the  process  being 

generated  is  equal  to  the  ratio  of  the  intensity  functions 

* * * 

evaluated  at  that  point,  i.e.  X(T^)/X  (T^) . 

More  formally: 

Theorem  1.  Consider  the  one-dimensional  non-homogeneous 

* * 

Poisson  process  {Nfc  : t ^ 0}  with  rate  function  X (t)  . 

* 

The  number  of  events,  Nt,,  in  the  fixed  interval  (0,t'] 

It  It 

has  a Poisson  distribution  with  parameter  y (0,t')  = u 
= A* ( t ' ) - A* (0) . 
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* * * * 

Let  T^,  Tj/  T^#  • Tn*  , be  the  times  of  the  events 

of  the  process  in  the  interval  (0,t']« 

Suppose  that  for  0 <_  t <_  t',  X (t)  <_  X*(t).  For 
* * 

i = 1,  2,  .../  Nfc,  delete  the  event  at  with  independent 
probability  1 - X (T*) /X* (T*) . 

Then  the  remaining  times  form  a non-homogeneous  Poisson 
process  with  rate  function  X(t)  in  the  interval  (0 , t ' ] . 

Proof : 

We  assume  that  X(t)  is  continuous  and  use  the  definition 
of  the  Poisson  process  based  on  incremental  probabilities. 
Thus  we  need  to  show  that  the  occurrence  of  an  event  in 
(t,t+dt]  is  independent  of  the  number  or  times  of  occurrence 
of  events  before  t,  and  that 

P{Nt+dt  " Nt  = 0}  = 1 " Mfc)dt  + o(dt), 

P{Nt+dt  " Nt  = = X(t)dt  + o(dt). 


and 


P{Nt+dt  - Nt  > 1}  = o(dt) . 
Now  we  have  that 
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P(no  event  from  {Nfc  : t ^ 0}  in  (t,t+dt]} 


★ 

= P{no  event  from  {Nfc  : t ^ 0}  in  (t,t+dt]}+  P{event 
* 

from  {N^  : t ^ 0}  in  (t,t+dt]  and  it  is  "thinned"} 

= 1 - A*(t)dt  + [A* (t) dt] • [1  - A (t) /A* (t) ] + o (dt) 

- 1 ” A* (t) dt  + A* (t) dt  - A*(t) -A(t)/A*(t) -dt  + o (dt) 

= 1 - A (t) dt  + o (dt) . 

Similarly: 

P{one  event  from  {Nfc  : t ^ 0}  in  (t,t+dt]  } 

•k 

= P{event  from  {Nfc  : t ^ 0}  in  (t,t+dt]  which  is  not  "thinned"} 
= A (t)dt*A(t)/A  (t)  + o(dt) 

= A ( t ) dt  + o ( t ) 

Also  it  follows  directly  that 

P(more  than  one  event  in  (t,t+dt]}  = o(dt) 

Moreover,  an  event  in  (t,t+dt]  is  independent  of  what 
happens  before  t because: 


i 


1.  (Nfc  : t > 0}  is  a Poisson  process  and  therefore  has 
independent  increments , and 

2.  The  thinning  uniform  random  variate  is  independent 

of  other  thinning  variates,  and  is  independent  of  the 

* 

Poisson  process  {Nfc  : t _>  0}. 

Q.E.D. 

Figure  1 shows  a graphical  representation  of  a particu- 
lar case  of  bounding  and  object  intensity  functions. 


i 
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Lng  Intensity  Function 


III.  ALGORITHMS  CONSIDERED 

A.  POISSON  DECOMPOSITION  AND  GAP-STATISTICS  ALGORITHM 
1.  Usage 

This  algorithm  is  the  one  found  by  PATROW  [Ref.  3] 
to  be  most  efficient  in  simulation  of  the  degree-two  exponen- 
tial polynomial  class  of  intensity  functions  and  its  imple- 
mentation by  PATROW  was  confined  to  that  group.  Basically, 

the  approach  is  to  decompose  the  intensity  function  (which 

2 

is  of  the  form  X(t)  = exp(a q + a^t  + a£t  ))  into  a degree- 
one  exponential  polynomial  function,  X^Ct)  = exp(bQ  + b^t) , 
and  a difference  function,  XD(t)  = X (t)  - X^ft).  The  points 
or  events  in  the  process  with  the  degree-one  exponential 
polynomial  function,  ( t ) , are  generated  over  the  interval 

(0,t' ] utilizing  the  computationally  fast  gap-statistics 
algorithm.  The  points  in  the  process  with  intensity  func- 
tion XD(t)  are  generated  using  conditioning-acceptance- 
rejection  techniques.  The  two  event  streams  are  then 
superposed  to  produce  the  event  stream  for  tne  non-homogeneous 
Poisson  process  with  the  intensity  function  X(t). 

In  the  case  where  X(t)  has  an  internal  maximum 
or  minimum  in  the  interval  (0,t']>  the  interval  is  par- 
titioned and  treated  as  two  separate  intervals  for  simula- 
tion with  the  event  streams  being  merged  in  the  final  step. 

Efficiency  is  optimized  by  maximizing  the  area 
under  the  degree-one  exponential  polynomial  intensity 
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function,  *L(t) , subject,  of  course,  to  the  constraint 
A^ft)  £ X(t)  for  0 <_  t <_  t'.  This  maximizes  the  use  of 
the  faster  gap-statistics  algorithm  and  minimizes  the  use 
of  the  conditioning-acceptance-rejection  algorithm  which 
is  slow  relative  to  the  gap-statistics  algorithm. 

PATROW  [Ref.  3]  deals  extensively  with  the  details 
of  this  algorithm  and  it  is  consequently  presented  here 
only  in  outline  form. 

2 . Algorithm  Statement 

a.  Categorize  the  intensity  function,  A (t) , into 

one  of  six  cases  by  examination  of  the  coefficients  a^ 

2 

and  a2  in  A(t)  = exp(aQ  + a-^t  + a2t  ).  Examples  of  each 
of  these  cases  are  shown  in  Figures  2 through  Figure  7 
located  in  Section  IV. 

b.  (1)  If  A(t)  is  monotone  increasing  or  monotone 
decreasing  over  the  interval  (Cases  I,  II,  IV  and  V;  see 
Figures  2,3,5  and  6),  decompose  A (t)  into  AT  (t) , which  is 
degree-one  exponential  polynomial,  and  Ap(t)  = A(t)  - AL(t). 
Thus  the  decomposed  functions  have  the  forms : 


AL(t)  = exp (bg  + bjt) 


Ad  ( t)  = exp  (Sq  + a1t  + a2t  ) - exp(bQ  + b^) 


AL(t)  + AD(t) . 


21 


Choose  k>0  and  so  as  to  maximize  the 
area  under  XL(t)  subject  to  XL(t)  <_  X(t)  for  all  t in 
(0 , t * ] . 

(2)  If  X (t)  is  not  monotone  over  the  interval 
(Cases  III  and  VI;  see  Figures  4 and  7) , partition  the 
interval  (0,t']  into  two  disjoint,  contiguous  subintervals, 
(0,b]  and  ( b , t ' ] - Choose  b as  the  (unique)  point  where 
X (b)  is  a maximum  (minimum)  of  X (t)  over  (0,t'].  Treat 
each  subinterval  as  in  b. (1) , applying  subsequent  steps  on 
each  subinterval  separately,  and  combining  results  as  the 
final  step. 

c.  Generate  points  in  the  Poisson  process  with 
degree-one  exponential  polynomial  intensity  function, 

XL(t),  using  the  gap-statistics  method. 

d.  Generate  and  order  points  in  the  Poisson  process 
with  intensity  function  XD(t)  using  the  conditioning- 
acceptance-rejection  method. 

e.  Merge  (superpose)  the  two  event  streams  from 
Step  3 and  Step  4.  The  merged  stream  is  from  the  non- 
homogeneous  Poisson  process  with  intensity  function  X(t). 

B.  THE  BASIC  THINNING  ALGORITHM 
1.  Usage 

The  thinning  theorem  is  implemented  in  a straight- 
forward manner.  Typically,  the  bounding  process  used  is 

* * 

homogeneous  Poisson  with  constant  rate  X , where  X is  an 
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upper  bound  of  X(t)  over  (0,t']«  In  this  case  efficiency 

* 

is  optimized  if  X is  the  least  upper  bound  (LUB)  of  X (t) 
over  ( 0 , t ' ] • 

2.  Algorithm  Statement 

a.  Generate  events  in  the  Poisson  process 

* * 

lNt  : t 0}  with  rate  function  X (t)  in  the  fixed  interval 

★ 

(0,t'].  If  the  number  of  events  generated,  n , is  such 

* 

that  n =0,  exit;  there  are  no  events  in  the  process 
{N  : t > 0}. 

* * * 

b.  Denote  the  (ordered)  events  by  T. , T_ , T *. 

12  n 

Set  i = 1 and  k = 0. 

c.  Generate  , uniformly  distributed  between  0 

and  1.  If  U.  < X(T*)/X  (T.  ) , set  k equal  to  k+1  and  T,  = T. . 

1 i>  X J\  1 

* 

d.  Set  i equal  to  i+1.  If  i <_  n , go  to  c. 

e.  Return  T^,  T2 , • Tn,  where  n = k,  and  also  n. 


C.  THE  ONE-AT-A-TIME  THINNING  ALGORITHM 
1 . Usage 

In  some  event-step  simulations,  it  is  customary  or 
necessary  to  generate  only  one  event  at  a time , rather  than 
an  array  of  events.  The  thinning  algorithm  is  easily 
modified  to  generate  the  next  event  in  the  non- homogeneous 
Poisson  process  with  intensity  function  X(t).  In  this  case, 
the  algorithm  utilizes  the  time  of  the  last  event,  T^_^; 
the  right  hand  limit,  t',  of  the  fixed  interval  over  which 

the  process  is  being  simulated;  and  the  bounding  process 

★ 

intensity  function,  X (t) . All  variates  are  generated 
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one  at  a time,  thus  no  arrays  are  required  for  storage. 


The  output  is  T^,  the  time  of  the  next  event,  if  any,  in 


the  interval  (T^_^,t']» 


The  algorithm  is  stated  here  for  the  case  in  which 

the  bounding  process  is  homogeneous  Poisson,  i.e., 

* * 

X (t)  = A , a constant  and  an  upper  bound  of  A(t). 
Specifically : 

T^  is  obtained  by  generating  and  cumulating 
★ ★ ★ 

exponential  (mean  = 1/A  ) random  variates  E.  ,,  E.  , , ...  , 

1 , x i , ^ 

for  i = 1,  2,  ...,  until  for  the  first  time. 


Di,j  x A(Ti-i  + Ei,i  + 


★ * 

+ Ei,j>A 


A detailed  algorithmic  statement  of  this  procedure 

follows . 

2.  Algorithm  Statement 

If  i = 1,  set  = 0 (i.e.  the  left  end  point 


of  the  interval) , otherwise,  T^_^  is  known.  Then  for  each 


i = 1 , 2 , ...  , the  time , T^ , of  the  event  in  the  non- 


homogeneous  Poisson  process  is  given  by  the  following: 

a.  Set  j = 1 


b.  Generate  E.  . , an  exponential  random  variate 
1 • 3 


with  mean  1/A 


* 2 * 

• ^ T.  i + l E.  j 

1 1 k=l 


is  greater  than  t'. 


exit;  there  are  no  more  points  in  the  interval 


c.  Generate  U ^ , uniformly  distributed  between  0 


and  1.  If 


i 


°jiA(Ti-l+  l Ei,k>/X' 

k*l 


set 


Ti  - Ti-i  + i Ei,k 

k=l 

and  exit. 

d.  Otherwise  set  j = j+1;  go  to  b. 

Note:  U.  . and  U.  are  uniformly  distributed 
1 / D J 

between  0 and  1. 
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IV.  METHODOLOGY  FOR  ALGORITHM  COMPARISON  f 

A.  MEASURES  OF  EFFECTIVENESS 

Two  quantifiable  measures  of  effectiveness  were  chosen 
as  yardsticks  for  algorithm  comparison.  These  were  compu- 
tational speed  and  computer  memory  requirements.  Some  other 
considerations,  such  as  programming  ease  and  robustness, 
are  discussed  in  Section  VI.  It  must  also  be  recalled  that 
the  classes  of  intensity  functions  for  which  the  two 
algorithms  are  usable  are  different.  The  Poisson  decomposi- 
tion and  gap-statistics  algorithm  is  only  easily  implemented 

for  a restricted  set  of  intensity  functions,  those  of  the 

2 

form  A (t)  = exp(aQ  + a^t  + a2t  ),  i.e.  the  degree-two 
exponential  polynomial.  Conversely,  the  thinning  algorithm 
is  valid  for  any  positive  intensity  function.  Thus  a direct 
comparison  can  only  be  made  in  that  subset  of  intensity 
functions  for  which  both  algorithm  implementations  are 
valid,  the  degree-two  exponential  polynomials. 

PATROW  [Ref.  3]  developed  six  sample  intensity  functions, 
all  special  cases  cf  the  degree-two  exponential  polynomial, 
and  these  are  used  herein  as  the  test  cases.  These  are 

described  in  Section  IV. C. 3 below.  \ 

1.  Computational  Speed 

Typically,  computer  time  is  a costly  commodity 
in  economic  terms.  It  may  also  be  a significant  factor 
in  determining  more  mundane  considerations  such  as  job 
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priority  and  thus  turn-around  time.  Thus  computer  run 
time  is  a natural  candidate  as  a measure  of  effectiveness 
for  comparing  competing  algorithms. 

PATROW  [Ref.  3]  utilized  a procedure  in  which 
event  streams  from  each  of  six  sample  intensity  functions 
were  replicated  several  times  in  "packages".  The  number 
of  replications  was  large  if  the  expected  number  of  events 
in  the  event  stream  was  small,  and  vice  versa.  Thus  the 
product  of  the  number  of  events  times  the  number  of  replica- 
tions was  kept  on  the  same  order  of  magnitude.  For  simplicity 
the  same  technique  was  used  here  although  results  showed  a 
wide  variation  in  the  run  times  for  the  six  packages. 

Programming  of  the  thinning  algorithm  was  done  so  as 
to  minimize  run  time  while  maintaining  parity  with  the 
Poisson  decomposition  and  gap-statistics  algorithm  wherever 
direct  comparison  could  be  made.  For  example,  shuffled 
random  numbers  are  called  in  both  programs  and  both  are 
dimensioned  to  accommodate  event  streams  of  up  to  5000  events. 

Undoubtedly  further  programming  refinements  exist 
which  might  increase  slightly  the  speed  of  one  or  the  other 
algorithm.  Also,  different  computers  might  have  unique 
features  which  could  be  exploited.  The  overall  purpose  here 
was  to  obtain  a relative  order  of  magnitude  comparison  and 
it  believed  that  this  objective  was  accomplished  in  every 
meaningful  sense. 


2.  Computer  Memory  Requirements 

This  is  the  second  obvious  means  of  comparing  two 
algorithms.  Again,  some  core  reduction  could  undoubtedly 
be  made  by  a sophisticated  programmer.  Most  notably,  core 
requirements  can  be  reduced  substantially  if  only  one  non- 
homogeneous  Poisson  variate  is  generated  at  a time  (the 
one-at-a-time  algorithm)  but  this  has  the  predictable 
effect  of  increasing  execution  time  considerably  (see 
Tables  IV,  V,  and  VI). 

B.  MEASUREMENT  CONSIDERATIONS 

Measurement  of  computer  memory  requirements  is  straight- 
forward and  deterministic. 

Measurement  of  computational  speed,  more  specifically 
Central  Processing  Unit  (CPU)  time,  is  quite  another  matter. 

First  of  all,  the  number  of  events  in  each  replication  of 
the  non-homogeneous  ^oisson  process  varies  causing  CPU  time 
to  be  a random  variable.  More  important,  however,  are  the 
effects  of  internal  computer  procedures. 

* 

In  the  first  place,  the  so-called  CPU  time  printed  out 
on  the  normal  IBM-360  output  has  only  a general  relationship 

I 'I  I 

to  the  actual  computational  time  required  by  the  CPU.  This 
is  caused  by  the  addition  of  certain  "overhead"  time.  This 
overhead  time  is  a function  of  the  number  of  other  programs 
in  the  system  as  well  as  such  factors  as  compilation  and 
linkage  times.  Thus  the  same  program  run  at  two  different 
times  may  differ  in  "CPU  time"  by  a factor  of  two  or  more. 

il 
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Program  execution  times  were  isolated  from  compilation 


and  linkage  times  by  the  use  of  a system  subroutine, 

GETIME.  This  subroutine  allows  the  user  to  initialize  an 
internal  timer  within  the  program  and  read  cumulative  time 
at  various  points  in  the  program.  Although  this  method  is 
not  exact,  it  does  measure  actual  elapsed  CPU  time  to  within 
a small  fraction  of  a second.  This  does  not,  however, 
entirely  alleviate  the  time-of-day  effect  experienced  when 
running  the  same  program  at  different  times.  That  is, 
although  the  elapsed  CPU  time  can  be  measured  accurately, 
the  same  program  will  generally  have  somewhat  different 
execution  times  each  time  it  is  run  [Ref.  4].  Theoretically, 
the  execution  times  would  be  constant  for  stand-alone  runs, 
i.e.  runs  with  no  other  competing  programs  in  the  system. 

This  is  rarely  realized  in  practice. 

These  considerations  lead  to  the  development  and  use  of 
the  side-by-side  setup  described  below.  This  method  appears 
to  be  statistically  sound  as  a means  of  dealing  with  the 
problems  of  time  measurement.  Due  to  the  differences  in 
execution  times  noted,  the  best  measure  of  effectiveness 
was  determined  to  be  a ratio  of  execution  times  for  the 
respective  algorithms,  rather  than  absolute  times. 

C . TEST  SETUP 

1.  Computational  Speed 

The  central  idea  here  was  to  equalize  the  effects 
of  non-essential  processes  on  each  algorithm.  This  was 
by  the  following  algorithm: 
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1.  Set  k = 1. 

2.  Zero  internal  time  clock. 

3.  Call  Algorithm  A.  Replicate  M times. 

4.  Read  internal  time  clock.  Store  time. 

5.  Zero  internal  time  clock. 

6.  Call  Algorithm  B.  Replicate  M times. 

7.  Read  internal  time  clock.  Store  time. 

8.  Set  k = k + 1.  If  k is  greater  than  k , go  to  9 . 

max 

Otherwise  go  to  2. 

9 . Compute  mean  and  variance  of  the  k execution  times 

max 

for  each  algorithm. 

10.  Compute  ratio  of  means. 

This  procedure  was  used  in  all  comparisons.  M, 

the  number  of  replications  per  package,  varied  between  30 

and  100  as  discussed  above.  K , the  number  of  times  each 

max 

package  was  replicated,  was  typically  set  equal  to  thirty. 
2.  Computer  Memory  Reauirements 


To  measure  computer  memory  requirements,  a small 
main  program  was  written,  calling  the  subroutine  which 
implemented  the  program  being  measured.  Total  program  length 
in  bytes  was  obtained  from  the  standard  computer  output  and 
the  core  alloted  to  the  main  program  was  subtracted  to 
obtain  the  desired  figure.  This  includes  all  library  rou- 
tines and  arrays  for  storage  of  event  times  and  arrays  of 
random  variates. 


The  core  requirements  are  deterministic  in  that 


they  do  not  change  from  one  run  to  another  but  are  strictly 
a function  of  the  program  coding. 

3 . Test  Cases  Utilized 

PATROW  [Ref.  3]  developed  six  sample  intensity 

functions  representing  the  possible  variations  in  sign  and 

relative  magnitude  of  the  coefficients  a^  and  a2  in  the 

2 

exponential  polynomial,  exp(ag  + a^t  + a2t  ). 

Since  the  sample  intensity  functions  were  designed 
to  test  different  aspects  of  the  Poisson  decomposition 
and  gap-statistics  algorithm,  they  were  also  used  for  com- 
parison here.  Although  each  algorithm  is  affected  by 
different  considerations,  the  test  cases  do,  coincidentally, 
put  the  thinning  algorithm  through  its  paces. 

For  continuity,  the  six  test  cases,  or  sample 
intensity  functions,  are  presented  below  in  Figures  2 through 
7. 
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SAMPLE  INTENSITY  FUNCTION 


Figure  7.  CASE  VI  - SAMPLE  INTENSITY  FUNCTION 


V.  EFFICIENCY  AND  PROGRAMMING  CONSIDERATIONS 


A.  GENERAL 

This  section  deals  with  factors  which  affect  the  per- 
formance of  the  thinning  algorithm.  Four  specific  areas 
are  presented  in  which  significant  gains  in  terms  of  com- 
putational speed  may  be  realized.  With  the  exception  of 
Section  V.D  (which  applies  only  to  the  case  of  exponential 
polynomial  intensity  functions)  these  considerations  apply 
to  the  general  class  of  intensity  functions. 

In  general  application,  one  of  the  primary  indicators 
of  efficiency  is  the  relative  size  of  the  area  under  the 
intensity  function  to  that  of  the  bounding  function,  i.e. 
the  ratio,  R,  given  by: 

t'  t* 

R = / A(t)dt //  A* ( t) dt 

0 0 

Since  both  numerator  and  denominator  are  simply  the  respec- 
tive integrated  rate  functions,  A(t),  evaluated  at  either 

end  of  the  interval,  R is  the  ratio  of  the  expected  number 

* 

of  events  in  the  two  processes,  i.e.  E [Nfc , ]/E [Nfc , ] . 

Case  1 of  the  sample  intensity  functions  is  particularly 
illustrative  (see  Figure  2) . The  intensity  function 
A(t)  = exp(1.6  + 0.015t  + 0.0005t2)  is  bounded  on  the 
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interval  (0,100]  by  a least  upper  bound  (LUB)  of  3294.47. 

If  a homogeneous  Poisson  process  with  rate  equal  to  the 

* 

LUB  is  used  as  the  bounding  function,  E [N^. , ] = 329,447 
points  will  be  generated  on  the  average.  Of  these,  all 
but  1464  will  be  rejected  on  the  average  (i.e.  E[Nt,]). 

The  ratio  of  the  respective  expected  values  is  thus 
1464/329,447  = 0.0044  = the  ratio  of  the  areas  under  the 
intensity  functions. 

Thus  a rough  relative  measure  of  the  efficiency  of  the 
thinning  process  in  a particular  situation  can  be  gained  by 
examining  a graph  of  the  two  intensity  functions,  even  if 


the  expected  values  are  not  easily  calculated.  This  pro- 
cedure may  also  be  an  indicator  in  deciding  whether  to 
partition  the  interval  and  use  different  bounding  functions 
on  each  subinterval. 

B.  UTILIZATION  OF  ARRAYS  OF  RANDOM  VARIATES 

Computer  generated  random  variates  are  used  both  in 

* 

generating  the  points  of  the  bounding  process  {Nfc  : t >_  0 } 
and  in  the  actual  thinning  process  itself.  Since  the 
number  of  variates  required  is  typically  large,  efficient 
generation  becomes  a programming  consideration  for  medium 
to  large  scale  simulations. 

The  basic  thinning  program  presented  herein  requires 
both  exponential  and  uniform  random  variates . Both  types 
are  obtained  utilizing  the  random  number  package,  LLRANDOM, 
developed  by  LEARMONTH  and  LEWIS  [Ref.  5].  Shuffled 
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tion.  Using  the  test  setup  of  Section  IV.  C,  average  times 

to  generate  varying  quantities  of  random  numbers  were 

determined.  Table  I reveals  the  relative  savings  realized 

by  calling  large  arrays  of  random  numbers.  Thus  considerable 

time  can  be  saved  by  generating  all  required  random  numbers 

from  one  subroutine  call.  Programming  difficulties  involve 

deciding  how  many  variates  to  generate.  The  general  goal 

is  to  generate  as  many  as  needed  while  keeping  the  unused 

excess  to  a minimum.  The  balance  used  was  to  generate  the 

expected  number  required  plus  an  excess  of  four  standard 

deviations.  For  example,  in  the  generation  of  the  bounding 

* * 

process,  the  expected  number  of  points,  E[Nt,]  is  a -t 

and  the  variance  is  the  same.  Thus  the  number  of  exponen- 

* ... 

tials  called  was  y + 4/y  where  y = X *t'.  Provision  is 

1 1 

. I ■ 

made  for  the  unlikely  (1  in  40,000)  case  that  more  are 

required . 

For  specific  applications  this  procedure  could  be 

improved  slightly.  For  example,  if  the  expected  number  of 
* 

points,  E[NfcI],  is  small,  e.g.  100,  then  the  expected 
excess  (four  standard  deviations)  comprises  forty  percent 
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of  the  total  whereas  for  large  E[Nt,],  e.g.  4000,  the 
expected  "waste"  is  only  about  six  percent.  In  the  former 
case,  reducing  the  "padding"  to  one  or  two  standard  devia- 
tions would,  on  the  average,  increase  efficiency  slightly 
although  the  probability  of  a second  subroutine  call  for 
more  random  variates  would  increase. 

As  an  example,  if  we  were  to  call  100  exponential 
variates  one  at  a time,  the  total  time  is  78,400  usee, 
compared  to  7,343  ysec  for  100  exponential  variates  called 
in  an  array.  For  100  + 4/100  = 140  variates,  the  time 
is  still  small  compared  to  78,400  ysec. 


Type  of 
Variate 

Number 

Called 

Total  Time 
(ysec) 

Mean  Time  Per 

Variate  (usee) 

Exponential 

1 

784 

784 

Exponential 

10 

1293 

129 

Exponential 

100 

7343 

73 

Exponential 

1000 

68046 

68 

Uniform 

1 

1213 

1213 

Uniform 

10 

1381 

138 

Uniform 

100 

3276 

33 

Uniform 

1000 

21544 

22 

Sample  Size  = 200  (each 

grouping) 

Table  I 

Generation 

Times  for 

Arrays  of  Shuffled  Random  Variates 

From  LLRANJOM 

C.  UTILIZATION  OF  INTENSITY  FUNCTION  LOWER  BOUND 


One  of  the  most  time  consuming  repetitive  operations 
is  the  computation  of  the  intensity  function  value,  A(t) , 
during  the  thinning  process.  In  the  case  of  the  exponential 
polynomial  intensity  function,  this  involves  one  power,  two 
multiplications,  two  additions,  and  one  exponentiation  for 
each  point  generated  in  the  bounding  process.  Since  points 
are  accepted  for  the  nortrhomogeneous  Poisson  process  when- 
ever the  uniform  (0,1)  chinning  variate  is  less  than  the 
★ 

ratio  A(t)/A  (t) , considerable  time  savings  result  if  the 

intensity  function  has  a positive  lower  bound,  say  A_, 

since  points  are  always  accepted  when  the  uniform  (0,1) 

* 

variate  is  less  than  the  ratio  A/A  . In  the  general  case, 
this  ratio  must  be  calculated  only  once.  The  expected  num- 
ber of  intensity  function  computations  which  are  alleviated 

* * 

by  the  use  of  the  lower  bound  is  given  by  (A/A  )E[Nfcl]  where 

★ 

Hs  a lower  bound  of  the  intensity  function;  A is  an  upper 

bound  of  the  intensity  function  (both  bounds  are  over  the 

* 

interval  (0,t'])  and  E[Nfcl]  is  the  average  number  of  points 
to  be  thinned,  i.e.  the  average  number  of  events  in  the 
bounding  process. 

It  is  clear  that  the  closer  A_  is  to  being  the  greatest 

* 

lower  bound  (GLB)  and  the  closer  A is  to  being  the  LUB, 
the  more  efficient  the  program. 

If  the  intensity  function  is  strictly  non-decreasing 
a further  (and  potentially  great)  improvement  is  realized 
by  initially  setting  A^  equal  to  A_,  and  then  setting  it 
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subsequently  equal  to  the  last  value  of  the  intensity 
function,  A(t).  This  results  in  a monotone  increasing 
lower  bound  and  thus  a decreasing  probability  of  evaluating 
the  intensity  function. 

Test  cases  II  through  VI  were  run  side  by  side  with  and 
without  the  use  of  a lower  bound  for  the  intensity  function. 

On  the  average,  the  program  which  did  not  utilize  a lower 
bound  required  twenty  percent  more  time  than  the  program 
using  a lower  bound.  Please  see  Appendix  B for  case-by- 
case comparison. 

D.  UTILIZATION  OF  EXPONENTIAL  VARIATES  FOR  THINNING  OF 
EXPONENTIAL  POLYNOMIAL  INTENSITY  FUNCTIONS 

The  time  requirements  for  evaluating  A(t)  were  discussed 

in  Section  V.C  above.  In  the  case  of  exponential  polynomial 

2 

intensity  functions,  e.g.  \(t)  = exp(aQ  + a^t  + a2t  ),  the 
major  contributor  to  computation  time  is  the  exponentiation 
operation.  Exponentiation  can  be  avoided  by  utilizing  the  } | 

following  relationship: 

* 

U.  \(t)/A  if  and  only  if 
E^  = -In  U^  _>  InA  - InA  (t)  = InA.  - (a^  + a^t  + a2t2)  , 
where : 

< 

U^  is  a uniform  (0,1)  random  variate, 

E^  is  an  exponential  random  variate  with  mean  one. 
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Thus  the  thinning  test  to  accept  points  from  the 
bounding  process  becomes: 

* 2 * 

If  >_  lnX  - (aQ  + a^t  + a2t  ) , for  t = T^ , accept 

* 

T^  as  a point  in  the  non-homogeneous  Poisson  process 

* 

with  rate  A(t);  otherwise,  reject  (i.e.  thin  it). 


The  key  to  this  relationship  lies  in  the  fact  that  if 
U is  distributed  uniform  (0,1),  then  -In  U is  distributed 
as  a unit  exponential  variate,  i.e.  an  exponential  variate 
with  mean  one.  This  is  shown  by  the  following: 

Let  U be  uniform  (0,1). 

Then  P{U  < x}  = P{ln  U < In  x)  s P{-ln  U >_  -In  x} 
but  P{U  _<  x}  = x,  thus  let  y = -In  x, 
then  P{-ln  U >_  y}  = exp(-y). 

Thus  -In  U is  distributed  as  a unit  exponential  variate. 

Although  more  time  is  required  to  generate  the  exponential 
random  variates  for  thinning  than  the  uniforms,  the  alle- 
viation of  the  exponentiation  operation  more  than  compen- 
sates for  the  additional  generation  time.  This  is  because 
SEXPON,  the  portion  of  LLRANDOM  which  generates  exponentials, 
generates  exponential  variates  by  the  Marsaglia  "rectangle- 
wedge-triangle"  method,  which  is  faster  than  taking  logarithms. 

Since  exponential  random  variates  are  used  in  the  genera- 
tion of  the  bounding  homogeneous  Poisson  process,  an  addi- 
tional time  savings  can  be  realized  by  using  the  variates 


ft  ^ 


which  are  left  over  (i.e.  not  used)  from  generating  the 
bounding  process  (these  are  generated  in  arrays) . 

For  the  test  cases  considered,  use  of  exponentials 
for  thinning  resulted  in  an  average  time  savings  of  ten 
percent.  Please  see  Appendix  B for  case-by-case  results. 

E.  RECYCLING  OF  THINNING  VARIATES 

As  mentioned  above,  a uniform  or  exponential  random 
variate  is  required  for  each  point  to  be  "thinned".  Each 
of  these  variates  requires  a significant  amount  of  time  for 
generation.  Obviously  a time  savings  would  be  realized  if 
fewer  variates  were  required. 

1.  Recycling  of  Uniform  Variates 

Assume  is  uniform  (0,1)  but  that  its  value  is 
unknown.  Assume  then  that  further  information  becomes 
available  that  is  less  than  a {0  <_  a <_  1)  , but  its  value 
is  still  unknown.  Then  is  uniformly  distributed  over 
the  interval  (0,a).  If  is  then  computed  by  "scaling 

up"  Ui,  i.e.  dividing  by  a,  then  Ui+1  is  uniform  (0,1). 
Similarly,  if  is  uniform  (0,1)  and  subsequent  information 
places  it  somewhere  above  a,  then  Ui+1  = (U^  - a)/(l  - a) 
is  uniform  (0,1).  Thus  by  conditioning  on  whether  the 
variate  is  greater  than  or  less  than  a given  value,  a new 
variate  can  be  computed  with  the  desired  properties. 

Moreover,  this  variate  is  independent  of  its  predecessor. 

In  the  thinning  algorithm,  each  point  is  tested 

* * 

using  a uniform  (0,1)  variate.  Specifically,  if  £ A(T^)/X 


the  point  is  accepted  as  a point  in  the  non-homogeneous 

* * 

Poisson  process.  Since  the  ratio  X(T^)/X  is  between  zero 

and  one,  and  the  only  test  is  whether  is  less  than  or 

* * 

greater  than  the  ratio  X(Ti)/X  , the  next  uniform  (0,1) 
variate,  can  be  generated  using  the  rules  above. 

The  algorithm  is: 

1.  Let  be  uniform  (0,1).  if  IL  is  less  than 
a = X (T*) /X* , let  U±+1  = U±/ (X (T*) /X*l ; exit. 

2.  Otherwise  let  U±+1  = [U±  - (X (T*)/X*]/ [1- (X (T*)/X*) ] 

Ui+1  is  uniform  (0,1). 

In  theory,  only  one  uniform  random  variate  is 

required  for  the  entire  thinning  process!  In  computational 

practice,  however,  care  must  be  exercised  because  of  the 

finite  capacity  of  the  computer  to  represent  numbers.  After 

ten  to  twenty  divisions  the  scaled  uniform  number  will 

consist  only  of  low-order  bits  of  the  random  number  and 

these  are  usually  not  uniformly  distributed. 

If  the  intensity  function  has  a positive  lower  bound, 

further  efficiencies  can  be  gained,  in  combination  with  the 

procedures  of  Section  V.C  above.  Since  multiplication  is 

* * 

computationally  faster  than  division,  the  value  1/(X/X  ) = X /X_ 

* 

can  be  precomputed  and  stored.  Thus  if  <_  X/X  , 

* 

Ui+1  = • X /\  can  be  computed  as  the  next  thinning  variate. 

Note  that  no  intensity  function  calculation  is  required. 

* 

However,  if  Ui  > \/\  , the  thinning  method  proceeds  with 

the  next  step,  evaluating  the  intensity  function  at  the 

* 

point  T^  and  determining  whether  or  not  to  thin  the  point. 
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Now,  further  information  is  known  about  U^.  Specifically, 
* * * 

either  > X(T^)/X  , in  which  case  TV  is  thinned,  or 

* ★ ★ 

X/X  < Ui  < X(T^)/X  . In  either  case  tV+1  can  be  computed 
by  "scaling  up"  IV. 


Thus,  the  algorithm  for  recycling  uniform  random 


variates  for  thinning  is  as  follows: 


1.  If  <_  X/\  , let  tV+^  = IV-  X /X_  and  exit. 

2.  If  X/X*  < ^ < X(T*)/X*,  let  Ui+1  = (X*-Ui-X>/(X(T*)-X) 


and  exit. 


3.  Otherwise,  U.  >X(T.  )/X  , let 


■ i" 

* * 


Ui+1  = {X  *ui-x<Ti) >/(*  - ^ (T±) ) 


By  precomputing  X /X_,  this  recycling  procedure 


requires  only  one  multiplication  in  the  case  where 


£ X/X  . Otherwise  one  multiplication,  two  subtractions 


and  one  division  are  required.  In  either  case  the  recycling 


procedure  is  generally  faster  than  generating  uniform 


random  variates  from  a random  number  generator,  even  when 


a logical  IF  statement  is  added  to  check  for  extreme 


values  ("small  bits") . 


2.  Recycling  Of  Exponential  Variates 


This  section  applies  only  where  the  intensity  func- 


tion is  exponential  polynomial.  Here  the  possibilities 


are  less  promising.  In  the  general  case  where  no  lower 


bound,  X_  is  used,  the  following  algorithm  would  apply: 


If  Ei  > In  X*  - In  X (T* ) , let  Ei+1  = Ei  - In  X*  + In  X(T*) 


Otherwise 


Let  Ei+1  = In  (X  - X(Ti))/(X  •exp(-Ei)  - X(T±)) 
where  E.  is  a unit  exponential  random  variate. 


~ _ . - w — 

★ it 

In  the  first  case,  E.  > In  a - In  X(T^),  a time 

★ 

savings  would  generally  be  realized  since  In  X could  be 

* 

computed  once  and  stored  and  In  X(T^)  is  simply  the  value 

* *2 

of  the  polynomial,  i.e.  aQ  + a1Ti  + a2Ti  , which  must  be 
computed  in  any  event.  In  the  second  case,  however,  the 
cure  is  truly  worse  than  the  illness.  It  is  faster  to 
simply  generate  another  exponential  variate,  assuming  they 
are  called  in  arrays. 

For  there  to  be  a time  savings,  however,  it  must 
be  possible  to  make  a reasonable  prediction  of  the  number 
of  exponentials  which  must  be  generated.  Otherwise  an 
excessive  number  of  calls  to  the  random  number  generation 
subroutine  may  destroy  the  gains  made  through  recycling. 

In  the  case  where  a lower  bound,  X_,  for  the  inten- 
sity function  exists  and  is  positive,  it  is  possible  to 
determine  the  expected  number  of  exponentials  which  must 
be  generated.  Variates  are  reused  if  they  are  greater  than 
ln(X*/X_).  That  is,  if  E^  > ln(X*/X),  then  Ei+1  = E^  - ln(X*/X^). 
Otherwise,  a new  (i.e.  non-recycled)  variate  is  used.  Thus 

it 

the  probability  of  not  recycling,  p,  is  X/X  and  the  number 

it  it 

of  variates  required  is  binomial  with  mean  n p,  where  n 
is  the  number  of  points  to  be  thinned. 

Empirical  results  for  the  five  test  cases  considered 
are  shown  in  Appendix  B.  Using  the  calling  rule  of  expected 
number  plus  four  standard  deviations  for  generating  thinning 
exponentials  yielded  inconclusive  results  as  compared  with 
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the  procedure  in  which  exponential  thinning  variates  are 
generated  in  arrays  with  no  recycling.  As  expected,  for 
larger  Nfcl  (Cases  II,  III  and  V),  recycling  provided  a slight 
time  advantage  (seven  percentage  maximum)  while  for  small 
Nfc , (Cases  IV  and  VI)  recycling  was  slower.  In  case  VI, 
recycling  caused  run  time  to  be  approximately  five  percen- 
tage greater  than  that  without  recycling.  Using  a calling 
rule  of  expected  number  plus  two  standard  deviations  reduced 
the  disadvantage  slightly  to  four  percent.  The  reason  that 
recycling  can  cause  longer  run  times  than  not  recycling  is 
that  an  additional  logical  IF  statement  is  required  for  the 
recycling  program.  Again,  when  exponentials  are  used  for 
thinning  and  the  mean  number  of  points  to  be  thinned  is 
on  the  order  of  two  or  three  hundred,  it  is  probably  not 
worth  the  effort  to  recycle.  If  several  thousand  "thinnings" 
are  required,  the  savings  may  indeed  be  worthwhile. 

Results  were  somewhat  surprising  in  the  general 
class  of  intensity  functions  for  which  uniform  variates  are 
used  for  thinning.  It  was  expected  that  a significant  savings 
would  be  realized  since  the  uniform  variates  can  be  recycled 
in  all  situations.  In  fact,  test  runs  were  run  in  which 
only  one  uniform  variate  was  generated  for  the  entire  run 
with  recycling  used  throughout.  The  only  program  statement 
which  added  time  was  a logical  IF  statement  to  preclude 
dividing  by  zero  and  to  diminish  the  probability  of  small 
bit  usage.  As  expected,  some  bias  was  experienced  in  the 
mean  and  an  unusually  high  variance  was  noted,  indicating  a 
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low  degree  of  "fidelity"  to  the  true  non-homogeneous  Poisson 
process  being  simulated.  Of  particular  interest  however, 
were  the  results  on  execution  time.  Since  this  setup 
essentially  gives  the  lower  bound  on  execution  time  for 
recycling  of  uniform  variates,  it  was  expected  that  signi- 
ficant time  savings  would  be  realized  in  comparison  to  no 
recycling.  In  practice,  however,  the  savings  were  minimal, 
with  a maximum  of  only  three  percent  savings.  This  is 
attributed  to  the  efficiency  of  the  LLRANDOM  package  in 
generating  uniform  variates  with  the  logical  IF  statement 
being  only  a secondary  cause  to  longer  execution  time. 

Appendix  B shows  results  for  both  exponential  and 
uniform  thinning  variates. 

The  key  point  to  keep  in  mind  here  is  that  the  above 
results  reflect  the  case  where  thinning  variates  are  called 
in  arrays,  i.e.  many  at  a time.  Thus  the  comparison  is 
between  a very  "fast"  variate  and  recycling.  As  discussed 
above,  calling  by  arrays  results  in  considerable  time  savings 
compared  to  one-at-a-time  generation  (up  to  fifty  times 
faster) . Thus,  in  the  case  where  a slower  random  number 
generator  is  utilized  or  where  variates  are  called  one  at 
a time,  use  of  the  lower  bound  may  indeed  result  in  con- 
siderable time  savings. 

F.  FINAL  PROGRAM 

A general  program  was  developed  which  incorporated  the 
efficiency  considerations  discussed  above.  The  program  is 


general  in  that  it  can  be  used  with  any  of  the  general 
class  of  intensity  functions,  whether  exponential  polynomial 
or  not.  The  program  is  essentially  four  programs,  each 
used  in  a specific  case.  The  program  classifies  simulations 
into  the  four  classes  by  asking  two  questions : 

1.  Is  the  intensity  function  exponential  polynomial? 

2.  Does  the  intensity  function  have  a positive  lower  bound? 

The  first  of  these  determines  whether  uniform  variates 

(general  case)  or  exponential  variates  (exponential  poly- 
nomial case)  are  used  for  thinning.  The  second  consideration 
merely  deletes  an  unnecessary  logical  IF  statement  in  the 
case  where  no  lower  bound  is  used. 

The  computer  program,  NHPP,  is  listed  after  the  appen- 
dices, and  requires  a user  supplied  subprogram  FUNCTION  FCN(T) 
to  compute  the  intensity  function  values  for  each  value  of 
t.  If  the  intensity  function  is  exponential  polynomial, 

only  the  exponent  portion  should  be  calculated,  i.e.  the 

2 

statement  FCN  = aQ  + a^t  + a2t  (for  degree-two  polynomial) 
should  appear  in  the  subprogram.  Otherwise  the  entire 
intensity  function  should  be  evaluated. 

Empirical  results  showed  that  the  final  program, 
utilizing  the  efficiency  considerations  mentioned  in  this 
chapter,  resulted  in  a program  which  ran  in  two-thirds  the 
time  of  the  basic  thinning  program. 

Please  see  Appendix  B for  case-by-case  results. 
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VI.  RESULTS,  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  GENERAL 

This  section  presents  the  results  of  comparison  of 
the  Poisson  decomposition  and  gap-statistics  algorithm 
with  three  variations  of  the  thinning  algorithm.  These 
variations  are  the  basic  thinning  algorithm,  the  modified 
thinning  algorithm  (final  program)  and  the  special  case 
one-at-a-time  algorithm.  Section  B presents  the  performance 
of  each  of  the  algorithms  when  measured  by  the  two  measures 
of  effectiveness,  computational  speed  and  computer  memory 
requirements.  Section  C examines  the  results  with  a view 
toward  identifying  the  strong  and  weak  points  of  each 
algorithm.  Section  D recommends  further  avenues  of  study. 

Again,  in  comparing  the  two  classes  of  algorithms,  one 
basic  distinction  must  be  kept  in  mind.  That  is  that  the 
Poisson  decomposition  and  gap-statistics  algorithm  as  imple- 
mented by  PATROW  [Ref.  3]  is  limited  to  a special  class  of 
intensity  functions,  i.e.  exponential  polynomial  of  degree 
two  (or  less) . Although  the  algorithm  could  be  adapted  to 
higher  order  polynomials  (by  further  bisection  of  inter- 
vals) , the  already  complex  programming  considerations  would 
grow  significantly.  In  contrast,  the  thinning  algorithm 
is  a completely  general  method  which  is  analytically  valid 
for  any  functional  form  of  permissible  intensity  functions 
(positive  and  right  continuous) . 
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The  results  presented  here  are  necessarily  limited  to 
that  class  of  intensity  functions  for  which  both  algorithms 
can  be  compared,  i.e.  the  degree- two  exponential  polynomial 
class.  The  purpose  was  to  determine  the  relative  performance 
of  the  thinning  algorithm  on  this  piece  of  common  turf  with 
the  heretofore  champion,  Poisson  decomposition  and  gap- 
statistics  . 

The  basic  result  is  that  the  thinning  algorithm  is 
indeed  quite  competitive  with  the  Poisson  decomposition  and 
gap-statistics  algorithm  in  the  area  of  mutual  validity. 

This,  combined  with  its  ease  of  programming  and  ability  to 
generate  variates  from  any  intensity  function,  make  the 
thinning  algorithm  a highly  attractive  tool  for  generating 
non- homogeneous  Poisson  processes  of  any  type. 

One  shortcoming  of  the  thinning  program  was  revealed  by 
the  first  test  case  considered  (see  Section  IV. C) . This  is 
a fast  rising  exponential  polynomial  which  rises  from  a 

value  of  five  to  almost  3300  over  the  interval  (0,100].  For 

★ 

A = 3294.47,  the  expected  number  of  points  in  the  bounding 
process  is  329,447  while  the  non-homogeneous  Poisson  process 
being  simulated  has  an  expected  number  of  only  1464  points. 

Thus  all  but  one  point  in  about  200  are  thinned  out.  The 
thinning  algorithm  could  be  more  efficiently  adapted  to 
this  case  by  partitioning  the  interval,  alleviating  the 
necessity  to  store  over  300,000  bounding  process  points. 

However,  the  efficiency  involved  would  still  be  low,  and  the 
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best  solution  appears  to  be  to  utilize  the  Poisson 
decomposition  and  gap-statistics  approach.  The  key  point 
here  is  that  the  problem  is  easily  recognized  beforehand, 
as  discussed  in  Section  V.A,  and  avoidable. 

Table  II  presents  a general  comparison  of  the  two 
algorithms . 


B.  MEASURES  OF  EFFECTIVENESS  RESULTS 

Chapter  four  details  the  comparison  procedure  utilized 
to  develop  the  following  results. 

1.  The  Basic  Thinning  Algorithm 

Salient  features  for  this  case  include  the  use  of 
uniform  variates  for  thinning  and  one-at-a-time  generation 
of  exponential  and  uniform  variates. 

Table  III  presents  the  results  for  each  of  the  five 
test  cases  run.  Algorithm  A is  computer  program  DEGTWO, 
the  Poisson  decomposition  and  gap-statistics  program 
developed  by  PATROW  and  listed  at  the  end  of  this  paper. 
Algorithm  B is  computer  program  NHPTHN , the  basic  thinning 
algorithm,  also  listed. 

The  thinning  algorithm  was  fastest  in  two  out  of 
the  five  test  cases  run  and  required  eighty  percent  of 
the  core  space  required  for  the  gap-statistics  algorithm. 

Table  VII  lists  core  storage  requirements  for  each 
algorithm. 

2 . The  Modified  Thinning  Algorithm  (Final  Program) 

This  section  compares  the  best  case  performance 

for  both  algorithms. 
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TABLE  III 

COMPARISON  OF  POISSON  DECOMPOSITION  AND  GAP- STATISTICS  AND  BASIC  THINNING  ALGORITHMS 


The  modified  thinning  algorithm  includes: 

1.  Use  of  exponentials  for  thinning  of  exponential 


polynomial  intensity  functions 

2.  Use  of  lower  bounds 

3.  Partial  recycling  of  exponential  thinning  variates 

4.  Use  of  exponential  variates  left  over  from  generation 
of  the  bounding  process. 

Each  of  these  refinements  are  discussed  in  Section  V. 

The  Poisson  decomposition  and  gap-statistics 
algorithm  used  was  again  the  implementation  by  PATROW 
[Ref.  3],  program  DEGTWO.  In  addition  to  the  normal  running 

of  the  program,  a second  set  of  comparisons  was  made  uti- 

* 

lizing  separately  calculated  values  for  c , the  bound  for 
the  conditioning-acceptance-rejection  routine.  This  is 
discussed  by  PATROW  [Ref.  3].  These  runs  are  indicated 
by  asterisk  (*)  in  Table  IV. 

In  the  first  case,  the  thinning  algorithm  was  faster 
in  four  out  of  the  five  test  cases,  with  the  best  relative 
performance  occurring  in  Case  VI. 

★ 

When  the  improved  values  of  c were  incorporated 
for  Cases  IV,  V,  and  VI,  the  Poisson  decomposition  and  gap- 
statistics  algorithm  improved  substantially,  winning  in 
three  out  of  the  five  test  cases. 

The  relative  advantage  in  computational  speed  was 
less  than  a factor  of  two  in  all  cases  with  a maximum  of 
1.83  to  1 in  favor  of  the  thinning  algorithm  in  case  VI. 
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Core  storage  for  the  thinning  algorithm  was  determined 
by  using  only  the  part  of  the  algorithm  which  used  exponen- 
tial thinning  variates.  This  precluded  the  requirement  for 
storing  5000  uniform  variates  which  are  not  used.  The 
Poisson  decomposition  and  gap-statistics  program  (DEGTWO) 
required  about  88,000  bytes  of  core  storage  as  compared  to 
about  94,000  for  the  thinning  program  (NHPP  modified  to 
exclude  unused  uniform  variate  array) . Detailed  results 
are  listed  in  Tables  IV  and  VI. 

3.  The  One-At-A-Time  Thinning  Algorithm 

As  discussed  in  Section  II. C,  the  one-at-a-time 

algorithm  was  developed  only  to  test  the  relative  efficiency 

of  the  algorithm  used  to  generate  the  next  event  in  a non- 

homogeneous  Poisson  process.  This  latter  requirement  may 

arise  in  event-step  simulation  where  only  the  next  event 

in  a non-homogeneous  Poisson  process  is  desired  rather  than 

an  array  containing  all  events  in  a specified  interval. 

Computationally,  the  one-at-a-time  algorithm  is 

quite  similar  to  the  basic  thinning  algorithm.  The  only 

essential  difference  is  that  the  basic  thinning  algorithm 

generates  and  stores  all  the  points  in  the  bounding  process 

★ 

(intensity  function  X ) before  thinning,  whereas  the  one- 
at-a-time  algorithm  generates  a point  in  the  bounding  process 
and  thins  that  point  before  continuing.  The  latter  method 
removes  the  requirement  for  an  additional  array  to  store 
the  bounding  process  points.  This  in  turn  saves  about 
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20,000  bytes  of  core  storage  requirement  when  the  programs 
are  dimensioned  to  accept  5000  points. 

As  implemented  here,  the  one-at-a-time  algorithm 
simply  generates  the  next  point  in  the  non-homogeneous  Poisson 
process  and  stores  it,  stopping  when  the  last  point  generated 
lies  outside  the  interval.  All  variates  in  this  program  are 
generated  one  at  a time.  The  results  shown  in  Table  V are 
thus  a good  indicator  of  the  relative  efficiency  of  using 
this  method.  As  can  be  seen,  the  one-at-a-time  algorithm 
(program  NHPOAT)  is  faster  than  the  Poisson  decomposition 
and  gap-statistics  algorithm  (program  DEGTWO)  in  three  of 
the  five  test  cases  run.  This  is  true  despite  the  fact  that 
DEGTWO  generates  all  variates  in  arrays,  taking  advantage 
of  the  time  economies  of  scale  mentioned  in  Section  V.  The 
one-at-a-time  algorithm  also  requires  forty  percent  less 
core  space. 

From  the  tables  one  can  also  see  that  since  both  the 
best  case  thinning  algorithm  (Table  IV)  and  the  one-at-a- 
time  thinning  algorithm  (Table  V)  are  compared  to  the 
Poisson  decomposition  and  gap-statistics  algorithm,  it  is 
possible  to  obtain  a reasonable  comparison  of  the  best  case 
thinning  algorithm  and  the  one-at-a-time  thinning  algorithm. 
For  example,  for  the  sample  intensity  function  used  in 
Case  II  (see  Figure  3),  the  ratio  (. 8127/. 5541)  = 1.47 
indicates  that  execution  time  for  the  one-at-a-time  thinning 
algorithm  is  almost  fifty  percent  greater  than  that  of  the 
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Case  E(N)  Algorithm  Replications  N o Time(Sec)  Ratio  of  Times 

Per  Package  Per  Package  (A  * B) 


COMPARISON  OF  POISSON  DECOMPOSITION  AND  GAP-STATISTICS  AND  ONE-AT-A-TIME  THINNING 

ALGORITHM 


! 


best  case  thinning  algorithm.  The  ratios  for  Cases  III, 

IV,  V and  VI  are  1.62,  1.27,  1.37,  and  1.21  respectively. 

The  tables  also  demonstrate  the  time-of-day  effect 
discussed  in  Section  IV.  That  is,  the  execution  time  for 
a given  program  is  not  the  same  each  time  it  is  run.  For 
example,  program  DEGTWO  took  21.8808  seconds  for  the  run 
recorded  in  Table  V compared  to  18.8917  seconds  for  the  run 
recorded  in  Table  IV.  For  this  reason,  ratios  of  execution 
times  were  chosen  as  the  measure  of  effectiveness  rather 
than  absolute  times. 

A FORTRAN  program  listing,  NHPNXT,  is  provided  at 
the  end  of  this  thesis.  This  program  generates  the  next 
point  in  a non-homogeneous  Poisson  process  with  a user 
supplied  intensity  subprogram,  FUNCTION  FCN.  This  program 
can  be  used  in  conjunction  with  event-step  simulation 
programs,  including  SIMSCRIPT,  where  it  is  desired  to  mini- 
mize core  space  (at  the  expense  of  speed)  or  where  only  one 
event  is  desired  at  a time.  Core  requirements  are  shown 
in  Table  VI. 

C.  CONCLUSIONS 

Both  the  thinning  algorithm  and  the  Poisson  decomposition 
and  gap-statistics  algorithm  include  two  general  types  of 
operations:  a "generating"  process  and  a "second  stage". 

For  the  Poisson  decomposition  and  gap-statistics  algorithm, 
the  generating  process  is  the  non-homogeneous  Poisson  process 
with  degree-one  exponential  polynomial  intensity  function. 

’ 

i 
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For  the  thinning  algorithm  (as  implemented  herein)  the 
generating  process  is  homogeneous  Poisson. 


I 


The  second  stage  for  the  Poisson  decomposition  and  gap- 
statistics  algorithm  is  the  actual  decomposition  and  genera- 
tion of  variates  by  the  conditioning-acceptance-rejection 
method.  For  the  thinning  algorithm,  the  second  stage  con- 
sists of  the  thinning  of  the  points  in  the  bounding  process. 
Thus  one  algorithm  generates  events  and  adds  more  events 
from  a second  process  while  the  other  generates  events  and 
subtracts  some  out. 

The  strongest  point  in  the  Poisson  decomposition  and 
gap-statistics  algorithm  is  the  highly  efficient  generation 
of  the  events  in  the  degree-one  exponential  polynomial  inten- 
sity function  process.  This  is  done  with  the  gap-statistics 
algorithm  which  is  two  to  five  times  faster  than  the  thinning 
algorithm  for  this  type  of  process  (see  Appendix  A) . At 
the  same  time,  the  conditioning-acceptance-rejection  routine 
is  relatively  quite  inefficient. 

There  are  many  considerations  in  predicting  the  relative 
success  of  the  two  algorithms,  i.e.  which  will  be  faster  in 
a given  situation.  For  example,  the  Poisson  decomposition 
and  gap-statistics  algorithm  is  affected  by  factors  such 
as  whether  or  not  partitioning  is  required,  the  percentage 
of  the  total  number  of  variates  which  come  from  the  degree- 
one  exponential  polynomial  process,  and  whether  time  rever- 
sal is  required.  For  the  thinning  algorithm,  the  fraction 
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of  the  lower  bound  divided  by  the  upper  bound  A/A  , would 
seem  to  be  a good  indicator  of  success. 

For  the  test  cases  considered,  however,  the  only  con- 
sistent indicator  was  the  expected  number  of  events  in  the  non- 
homogeneous  Poisson  process  being  simulated.  The  smaller  the  number 
of  events  in  the  non-homogeneous  Poisson  process  being  simu- 
lated, the  better  the  relative  performance  of  the  thinning 
algorithm  over  the  Poisson  decomposition  and  gap-statistics 
algorithm.  Thus  it  appears  that  each  algorithm  has  a fixed 
and  variable  part  in  terms  of  time.  The  thinning  algorithm 
has  a shorter  "setup"  cost  in  terms  of  time  but  the  variable 
cost  or  "cost  per  additional  variate"  seems  to  be  smaller 
for  the  Poisson  decomposition  and  gap-statistics  algorithm. 

The  exact  cause  of  this  phenomenon  is  not  known  although 
it  appears  to  be  centered  in  the  conditioning-acceptance- 
rejection  routine. 

In  the  larger  spectrum  of  non-homogeneous  Poisson  process 
generation,  it  seems  clear  that  the  thinning  algorithm  is 
the  best  all-around  method  available. 


D.  RECOMMENDATIONS 

Two  specific  areas  for  further  study  are  recommended. 

First,  the  thinning  algorithm  as  implemented  here  uses 
only  homogeneous  Poisson  processes  for  bounding.  It  might 
be  worthwhile  to  investigate  the  possibility  of  using  other 
processes,  such  as  non- homogeneous  Poisson  processes  with 
degree-one  exponential  polynomial  intensity  functions,  as 


bounding  processes.  This  would  allow  the  efficient  gap- 
statistics  algorithm  to  be  utilized  although  function 
evaluation  would  in  general  become  more  time  consuming. 

The  second  area  is  in  finding  the  optimum  method  for 
generating  the  degree-two  exponential  polynomial  class  of 
intensity  function.  These  will  undoubtedly  remain  of 
interest  due  to  their  statistical  properties.  Here,  it 
seems  clear  that  the  best  features  of  the  two  algorithms 
can  be  combined.  Specifically,  the  Poisson  decomposition 
and  gap-statistics  algorithm  can  be  modified  to  use  thinning 
rather  than  conditioning-acceptance-rejection  for  generating 
the  points  in  the  difference  function  process,  i.e.  the 
process  with  intensity  function  X-  - X (t)  - X_  (t)  . Also 
the  criterion  for  the  decomposition  might  preferably  be 
that  the  intensity  function  of  the  remainder  be  monotonically 
increasing.  This  would  make  it  easy  to  find  the  upper 
bound  for  the  function,  and  the  most  efficient  version  of 
the  thinning  algorithm,  that  where  the  number  of  computations 
of  the  intensity  function  is  minimized,  could  be  used. 
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APPENDIX  A 


r 

■ 

GENERATION  OF  DEGREE-ONE  EXPONENTIAL 
POLYNOMIAL  INTENSITY  FUNCTIONS 

The  generating  process  for  the  Poisson  decomposition 
and  gap-statistics  algorithm  is  a non- homogeneous  Poisson 
process  with  degree-one  exponential  polynomial  intensity 
function.  This  is  generated  by  using  the  gap-statistics 
method,  which  is  subroutine  NHPP2  in  the  DEGTWO  program 
(see  listing  below) . 

II 

To  determine  the  relative  speed  of  the  thinning  algorithm 
compared  to  the  gap-statistics  algorithm,  two  simple  degree- 
one  exponential  polynomial  intensity  functions  were  developed 
and  simulated. 

Table  VII  presents  the  results.  Case  A is  a monotone 
decreasing  intensity  function,  X (t)  = exp (3. 4 - 0.02*t) 
over  the  interval  (0,100].  Case  B is  a monotone  increasing 
intensity  function,  A(t)  = exp (.693  + 0.03-t)  over  the 
interval  (0,50]. 


Results  show  that  the  gap-statistics  algorithm  is  from 
two  and  a half  (Case  B)  to  four  and  a half  (Case  A)  times 
faster  than  the  thinning  algorithm. 
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APPENDIX  B 

RESULTS  OF  EFFICIENCY  MODIFICATIONS 

This  appendix  presents,  in  tabular  form,  the  results 
of  comparison  of  the  programming  modifications  listed  in 
Section  V. 

Table  VIII  shows  the  effects  of  utilizing  lower  bounds 
for  the  intensity  function.  Test  conditions  include: 

1.  Use  of  uniform  thinning  variates 

2.  Recycling  of  thinning  variates 

3.  Use  of  arrays  of  variates 

Table  IX  shows  the  gains  realized  by  employing  exponen- 
tial thinning  variates  in  contrast  to  uniform  thinning 
variates  for  exponential  polynomial  intensity  functions. 
Test  conditions  include: 

1.  Use  of  arrays  of  variates 

2.  Recycling  of  thinning  variates 

3.  Use  of  lower  bounds  for  intensity  function 

Table  X shows  the  results  of  recycling  versus  no 

recycling  where  uniform  variates  are  used  for  thinning 
while  Table  XI  shows  the  same  comparison  when  exponential 
variates  are  used  for  thinning.  For  both  cases,  test 
conditions  include: 

1.  Use  of  arrays  of  random  variates 

2.  Use  of  lower  bounds  for  intensity  function 
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Case  I E(N)  Algorithm  I Replications  M*  N Time(sec)  I Ratio  of  Times 


LOWER  BOUND  VERSUS  NO  LOWER  BOUND;  UNIFORM  THINNING  VARIATES 


Case  E(N)  Algorithm  Replications  M*  N Time(sec)  Ratio  of  Times 

Per  Package  Per  Package  (At  B) 


EXPONENTIAL  THINNING  VARIATES  VERSUS  UNIFORM  THINNING  VARIATES  FOR 
EXPONENTIAL  POLYNOMIAL  INTENTISY  FUNCTION 


Case  E(N)  Algorithm  Replications  M*  N Time(sec)  Ratio  of  Times 

Per  Package  Per  Package  (A  t B) 
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RECYCLING  VERSUS  NO  RECYCLING;  EXPONENTIAL  THINNING  VARIATES 


Table  XII  presents  the  results  of  incorporating  all 
of  the  programming  improvements  into  program  NHPP.  The 
final  thinning  program,  NHPP,  is  compared  to  the  basic 


thinning  program  without  modifications,  NHPTHN . The 
essential  differences  are: 


NHPP  (final  program) 
Arrays  of  variates  generated 


Exponential  variates  used 
for  thinning 

Lower  bound  of  intensity 
functions  utilized 

Thinning  variates  recycled 


NHPTHN  (basic  program) 

Variates  generated  one 
at  a time 

Uniform  variates  used 
for  thinning 

Lower  bound  = 0.0 


No  recycling  used 
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FINAL  THINNING  PROGRAM  VERSUS  BASIC  THINNING  PROGRAM 


ooooooooooooooooooooooooooooooooooooonoooooooooooooooooooooooooonnooooo 


THIS  PACE  IS  BEST  QUALITY  PRACTICABLE 
ZBOtf  OOPY  7UMUSHED  TO  DOC 


SUBROUTINE  NHPP 
SUBROUTINE  NHPP 
PURPOSE 

SIMULATES  A NON- HOMOGENEOUS  POISSON  PROCESS 
WITH  INTENSITY  FUNCTION  FCN(XI  OVER  A GIVEN 
INTERVAL  USING  THE  THINNING  ALGORITHM. 


USAGE 

CALL 


NHPP(IS»EL»ER»UB»XMIN»NTYPE»N,IER) 


DESCRIPTION  OF  PARAMETERS 

IS  ' - 


ANY  INTEGER  WITH 


EL 

ER 

UB 


X MIN 


RANDOM  NUMBER  SEED. 

NINE  OR  LESS  DIGITS. 

LEFT  END  POINT  OF  INTERVAL 
RIGHT  END  POINT  OF  INTERVAL 
UPPER  BCUNO  OF  THE  INTENSITY  FUNCTION 
FCN(X)  OVER  THE  INTERVAL  (EL  » ER)  • THE 
CLOSER  UB  IS  TO  THE  LUB,  THE  MORE 
EFFICIENT  THE  PROGRAM 

- LOWER  BOUND  OF  THE  INTENSITY  FUNCTION 
OVER  THE  INTERVAL.  THE  CLOSER  XMIN  IS  TO 
THE  GLB  , THE  MORE  EFFICIENT  THE  PROGRAM. 
NTYPE  - 1 IF  THE  INTENSITY  FUNCTION  IS 

EXPONENTIAL  POLYNOMIAL.  I.E.  OF 
THE  FORM  EXP (A  + A *X  ♦ A2*X**2  +...) 
0 OTHERWISE 

N - THE  TOTAL  NUMBER  OF  EVENTS  IN  THE 
NON-HOMOGENEOUS  POISSON  PROCESS 
IER  - ERROR  FLAG.  IER  HAS  FOLLOWING  MEANINGS: 
1.  ..ER  IS  LESS  THAN  EL 

2.. .UB  IS  NON-POSITIVE 

3. . .XMIN  IS  NEGATIVE 

4..  .MORE  THAN  5000  EVENTS  REQUIRED  FOR 

BOUNDING  PROCESS:  STORAGE  CAPACITY 
EXCEEDED. 

5..  .XMIN  IS  GREATER  THAN  UB 


COMMENTS 

CALLING  PROGRAM  MUST  HAVE  A COMMON  REGION, 
DCNNA,  OF  DIMENSIQN(5000) 

EXAMPLE:  DIMENSION  T( 5000) 

C3MM0N/D0NNA/T 

TIMES  TO  EVENTS  WILL  BE  STORED  IN 
CELLS  T(l)  THROUGH  TIN) 

THE  INTENSITY  FUNCTION  IS  USER  SUPPLIED. 

IF  THE  INTENSITY  FUNCTION  IS  NOT  EXPONENTIAL 
POLYNOMIAL,  I.E.  IF  NT YP  E * 0,  THE  ENTIRE 
INTENSITY  FUNCTION  IS  EVALUATED. 

EXAMPLE:  FUNCTION  FCN(X) 

A = 1.0 
FCN  = B *X 
RETURN 
END 

IF  THE  INTENSITY  FUNCTION  IS  EXPONENTIAL 
POLYNOMIAL,  I.E.  NTYPE  * 1,  ONLY  THE 
POLYNOMIAL  IS  EVALUATED. 

EXAMPLE:  FUNCTION  FCN(X) 

A a 1 .0 
A1  * 0.5 
A2  3 0.05 

FCN  = A ♦ A1*X  * A2*X**2 
RETURN 


! ' 
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00  5 I S1 1 NS  T A R «v.«ffYnBaiSHH)IO 

IF  (EB.LT.BND)  GO  TO  2 JSOI 00PY  TURAISHHu  iv 

K * K+l 

TT(K»  = TIMESd) 

EB  * EB-BNO 
GO  TO  5 

2 VA L = -FCN  (TI  MES  ( I ) )+UBLN 
IF  (EB.LT.VAL)  GO  TO  3 

K « K+l 

TT(  K)  = TIMESd  > 

3 KK  * KK  + 1 

CHECK  TO  SEE  IF  MORE  UNIT  EXPON  NEEOEO 

IF  (Kk.LE.NCHK)  GO  TO  4 

GENERATE  MORE  UNIT  EXPONENTIALS  FOR  THINNING 

KKK  = KKK+KK 

FN  = P*FL  OAT ( NST  AR-KKK ) 

NEB  = MAX0(1,IFIX( PN+4.0*SQRT(PN*CI) > 

CALL  SEXPON  <IS,EE,NE3) 

KK  a 1 
NCHK  = NEB 

4 EB  = EECKK) 

5 CONTINUE 
N * K 

GO  TO  14 

LOG  QUADRATIC  WITH  NO  MINIMUM 

6 CONTINUE 

CALL  REORD  (N  STAR.N  EXP  ,NLEF  T) 

NCALL  = NSTAR-NLEFT 
IF  (NCALL. LE.O)  GO  TO  7 
CALL  SEXPON  llSt  EE,  NCALL) 

SET  VARIABLE  S 

7 K * 0 

UBLN  = ALOG(UB) 

I COUNTS  HPP  EVENTS 
K COUNTS  NHPP  EVENTS 

00  e I =1 1 NSTAR 
VAL  = -FC  N (TI  MES  ( I ) ) +UBLN 
IF  (£E( I ) .LT.VAL)  GC  TO  g 
K = K+l 

TT  (K  ) = TIMESd) 

8 CONTINUE 
N = K 

GO  TC  14 

INTENSITY  FUNCTION  IS  NOT  LOG  QUADRATIC 

DOES  IT  HAVE  A MINIMUM  CR  IS  MINIMUM  LESS  THAN 
PCTMN  Of  MAX? 

9 PCT  = XMIN/UB 

IF  (PCT  «i.T  . PCTMIN  ) GO  TO  12 

USE  MIN  IMUM 

INITIALIZE  VARIABLES 

NUMF  * NSTAR 

CALL  SRAND  ( IS»U, NUNIF ) 

K « 0 

FF  « 1.0/UB 
I COUNTS  HPP  EVENTS 
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c 
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10 


11 


DO  11  I*1.NSTAR 
IF  (U(I). GT.PCT)  GO  TO  10 
K = K+l 

TT(KI  * TIHESd  ) 

GO  TO  11 

VAL  - FCNiTiNESII  ) ) *FF 
IF  (Utll.GT.VALl  GO  TO  11 
K * K+i 

TTCKJ  * T 1 M£S (II 
GO  TO  11 
CONTINUE 
N * K 
GO  TO  14 


NO 

12 


MINIMUM 


NUNIF  « NSTAft 

CALL  SRAND  ( IS, U, NUNIF) 

K * 0 

FF  * 1.0/UB 

DO  13  1*1 ,N  STAR 

VAL  * FC  N(T IMES ( I ) )*FF 

IF  (UllJ.GT.VAL)  GO  TO  13 


ACCEPT  POINT 


K * K+l 

TTiKI  = TIMES  (I  ) 


REJECT  PGINT 


13 


CONTINUE 
N * K 
14  RETURN 
END 


SUBROUTINE  HPP 


SUBROUTINE  HPP  GENERATES  PCINTS  IN  A HOMOGENEOUS 
POISSON  PROCESS  WITH  INTENSITY  FUNCTION  * UB 


SUBROUTINE  HPPU  S , EL,ER,UB,  NTYPE ,NST AR ,NL E FT , N EXP, I ER  ) 
DIMENSION  T I M ES  ( 5000) , EE15000) 

COMMON  /ANNE/  TIMES, EE 


INITIALIZE  VARIABLES 


EXMEAN  * 1.0/LB 

NEXP  * IF  IX  (U 8* ( ER  - EL  ) +4 .0 *SQRT  (U8*(  ER— EL  ) ) ) 
IF  (NEXP.  GT.  5000)  NEXP=5000 
CALL  SEXPON  (IS, EE, NEXP) 

SUM  = EL 
ISTART  = 0 
ISTOP  * NEXP 


DO  1 JJ*1 ,NEXP 
E * EXMEAN*EE  (JJ  ) 

SUM  * SUM+E 

TIMES(JJ)  * SLM 

IF  (SUM.LT.ER  ) GO  TO  1 

NSTAR  * JJ-1 

C-0  TO  5 

CONTINUE 


EXCEPTIONAL  NUMBER  OF  POINTS  NEEDED 


2 ISTART  * I START+J J 

NEXP  * IFIX(UB*(ER-SUM  )«-4.Q*SQRT(  UB*(ER-SUM  ) ) ) 
IF  (NEXP.EQVOI  NEXP-1 
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JSTCP  = I ST ART+NEXP 
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GO  TO 
* NEXP ) 


IF  (IST0P.GT.5000) 

CALL  SEXPON  (IStEE 
CO  3 J J*I » N EXP 
E * EXMEAN*EE (JJ) 

SUM  = SUM+E 

KK  * JJ+I  START 

TIMES(KK)  * SUM 

IF  (SUM.LT.ER)  GO  TO  3 

NSTAR  = KK— I 

GO  TO  5 

CONTINUE 


GO  TO  2 


MORE  THAN  5000  EXPONENTIALS  NEEDED 


4 IER  « 4 
GO  TO  6 


CALCULATE  NUMBER  OF  EXPONENTIALS  NOT  USED 


5 NLEFT  i 

6 RETURN 
ENC 


NEXP-JJ 


• ..SUBROUTINE  REORD 

SUBROUTINE  REORD  REORDERS  EXPONENTIALS  LEFT  OVER  FROM 
SUBROUTINE  HPP  FOR  USE  AS  THINNING  VARIATES  IN  NHPP 


SUBROUTINE  RECRD  ( NCHK  ♦ N EXP  ,NL  EFT  ) 
DIMENSION  EE1  5000J.  TIMES(5000) 
COMMON  /ANNE/  TIMES, EE 


ARE  ENOUGH  EXPONENTIALS  LEFT  OVER  FRCM  HPP? 

IF  (NLEFT. GE. NCHK)  GO  TO  4 
MORE  EXPONENTIAL  S NEEDED 
REORDER  TOP  TO  BOTTOM  OR  VICE  VERSA? 

IF  (NEXP.  LT. NCHK)  GO  TC  2 
RECROER  BOTTOM  TO  TOP 


N1  = NC  HK  -NLEFT 
N2  a NEXP-NLEFT 


CO  I 1=1,  NLEFT 

J = Nl+I 

K = N2+I 

EE ( J ) = EE(K) 

CONTINUE 

GO  TO  6 


REORDER  TOP  TO  BOTTOM 


NCHKH 

NEXP+1 


8 

C 


NCHKO 
NE  XPO 
DO  3 1=1  .NLEFT 
J = NCHKO— I 
K = NEXPO-I 
EE(J)  * EE(K) 
CONTI  NUE 
GO  TO  6 


ENOUGH  LEFT  OVER  FROM  HPP  FOR  ALL  THINNING 
4 


N1  = NEXP-NLEFT 
DO  5 1=1, NCHK 
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J * Nl+I 

EE  (I)  * E E ( J 1 

CONTINUE 

RETURN 

END 
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SUBROUT INE  NHPTHN 


PURPOSE 

SIMULATES  A NONHOMOGENEOUS  POISSON  PROCESS 
WITH  INTENSITY  FUNCTION  FCN(X)  USING  . 

THE  THINNING  ALGORITHM 

USAGE 

CALL  NHPTHN(IS,EL,ER,UB,N,IER  1 


DESCRIPTION  OF  PARAMETERS 


IS  - RANDOM  NUMBER  SEED.  ANY  INTEGER  WITH  NINE 
OR  LESS  OIGITS. 

EL  - LEFT  END  POINT  OF  INTERVAL 

ER  - RIGHT  END  POINT  OF  THE  INTERVAL 

UB  - UPPER  BOUND  OF  THE  INTENSITY  FUNCTION  , FCM X ) » 
OVER  THE  INTERVAL  (EL,ER).  THE  CLOSER  UB  IS 
TO  THE  LEAST  UPPER  BOUND, LUB,  THE  MORE 
EFFICIENT  THE  PROGRAM.  UB  MUST  BE  STRICTLY 
POSITIVE 

N - THE  TOTAL  NUMBER  OF  EVENTS  IN  THE 
NQN-HOMOGENEOUS  PCISSON  PROCESS. 

IER  - ERROR  FLAG.  IER  HAS  FOLLOWING  MEANINGS*. 

1.  ..ER  IS  LESS  THAN  EL 
2...UB  IS  NON-POSITIVE 


COMMENTS 

CALLING  PROGRAM  MUST  HAVE 
SCOTT,  OF  DI  ME  NS  I ON  (5000) 


COMMON  REGION, 


EXAMPLE*.  DIMENSION  T (5000  ) 

C GMMON / SCOTT/T 

TIMES  TO  EVENTS  WILL  BE  STORED  IN  CELLS 
T(  I ) THROUGH  TIM). 


PROGRAMMER: 


LCDR  JOHN  SCOTT  REDD,  USN 
AUGUST  1573 


SUBROUTINE  NHPTHN  (I S , EL, ER,U B , N , I ER ) 

DIMENSION  TIMES  ( 5000)  , TTTOOOO) 

COMHCN  /SCOTT/  TTT 
EXTERNAL  FCN 
CALL  OVFLOW 
INITIALIZE  VARIABLES 
IER  = 0 

IF  (EL.GE.ER)  IER  = 1 
IF  (UB.LE.O.O)  IER  = 2 
Ip  (IER.NE.O)  RETURN 

GENERATE  POINTS  IN  HOMOGENEOUS  PCISSCN  PROCESS  WITH  RATE 
= U8  . 

EXMEAN  = 1.0/LB 
1 = 1 
SUM  = EL 

1 CONTINUE 

CALL  EXPON  ( I S , E,  1) 

E = E+EXMEAN 
SUM  =*  SUM+E 

IF  (SUM.GT.ER)  GO  TO  2 
TI MES(I)  = SUM 
I = I+l 
GO  TO  1 

2 CONTINUE 
NSTAR  * 1-1 

COMMENCE  THINNING  THE  NSTAR  POINTS 
FI  « 1.0/UB 
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K * 0 

IF  (NSTAR.EQ.QJ  GO  TQ  4 

00  3 I*1,MSTAR 
CALL  SRAND  (IStUtl  I 
RATIO  » FCN(TIMES(I))*F1 
IF  (U.GT.RATIO)  GO  TO  3 
K » K+l 

„ TTT  t K ) * TIMES(I) 

3 CONTINUE 

N = K 
RETURN 

4 N * C 
RETURN 
ENC 
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SUBROUTINE  NHPOAT 


SUEROUHNE  NHPOAT 


am  MSB  IS  bust  wilt"  moncisu 
root  oa ri  fukiushed  to  ddc  — 


purpose 

SIMULATES  A NON-HOMOGENEOUS  POISSON  PROCESS 
WITH  INTENSITY  FUNCTION  FCN(X>  OVER  THE 
INTERVAL  ( EL  »ER  ) USING  THE  ONE-A  T-A-TI  ME 
THINNING  ALGORITHM 

...usage 

CALL  NHPOAT(IS,EL,ER,UB,N,IER» 


...DESCRIPTION  OF  PARAMETERS 


IS  - RANDOM  NUMBER  SEED.  ANY  INTEGER  WITH  NINE 
CR  LESS  DIGITS . 

EL  - LEFT  END  POINT  OF  INTERVAL 
ER  - RIGHT  END  POINT  OF  THE  INTERVAL 
UB  - UPPER  BOUND  OF  THE  INTENSITY  FUNCT  ION,FCN(X), 
OVER.  THE  INTERVAL  (EL, ER  ) • THE  CLOSER  UB  IS 
TO  THE  LEAST  UPPER  BOUND,  LUB  , THE  MORE 
EFFICIENT  THE  FROG  PAM.  UB  MUST  BE  STRICTLY 
POSI  TI  VE 

N - THE  TOTAL  NUMBER  OF  EVENTS  IN  THE 
NON-HOMQGSNEOUS  PCISSON  PROCESS. 

IER  - ERROR  FLAG.  IER  HAS  FOLLOWING  MEANINGS; 

1.. .ER  IS  LESS  THAN  EL 

2.. .UB  IS  NON-POSITIVE 
...COMMENTS 

CALLING  PROGRAM  MUST  HAVE  A COMMON  REGION,  DONNA,  OF 
DIMENSION! 5C00) 

EXAMPLE;  DIMENSION  TTI5000I 
COMMON/DONNA /TT 

TIMES  TO  EVENTS  WILL  BE  STORED  IN  CELLS 
T ( 1 ) THROUGH  T(N) 

SUBROUTINE  NHPOAT  (IS,EL,E'R,U3,N,IER) 

DIMENSION  TT(SOOO) 

COMMON  /DONNA/  TT 
EXTERNAL  FCN 
CALL  CVFLOW 

...INITIALIZE  VARIABLES 


IER  * C 
IF  (EL.GE.ER)  IER  = 1 
IF  (UB.LE.O.O)  IER  * 2 
IF  (IER.NE.Ol  RETURN 


1*1 
EXMEAN  * 
SUM  * EL 


l.O/UB 


GENERATE  POINT  IN  BOUNDING  PROCESS 

1 CONTINUE 

CALL  EXPON  ( IS  , E,  1 » 

E * E*EXMEAN 
SUM  = SUM+E 

IF  <SUM.GT.ER)  GO  TO  2 

THIN  THE  POINT 

CALL  RANDOM  (IS,U,1) 

RATIO*  FCN  ( SUM  ) *EXM  EAN 
IF  (U.GT.RATIC)  GO  TO  1 
TT  ( I ) * SUM 
I * I+l 
GO  TO  1 

2 N * 1-1 
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SLBROUT 1NE  NHPNXT 


SUBROUTINE  NHPNXT  GENERATES  THE  TIME  OF  THE  NEXT 
EVENT  IN  A NGN-HOMOGENEOUS  POISSON  PROCESS  WITH 
RATE  FUNCTION  FCN(X)  (USER  SUPPLIED). 


USAGE : 

CALL  NHP NXT( I S, UB, GLB, XLAST, ER, XNEXT, IER) 


DESCRIPTION  OF  PARAMETERS? 


IS 


UB 


GL  E 


XLAST 

ER 

XN  EXT 


I ER 


RANDOM  NUMBER  SEED.  ANY  INTEGER  WITH 
NINE  OR  LESS  DIGITS. 

UPPER  BOUND  OF  THE  INTENSITY  FUNCTION 
OVER  THE  INTERVAL  (XLAST.ER). 

GREATEST  LOWER  BOUND  OF  THE  INTENSITY 
FUNCTION  OVER  THE  INTERVAL  (XLAST.ER). 

SET  = 0 IF  UNKNOWN. 

THE  TIME  OF  THE  LAST  EVENT  IN  THE  PROCESS 
RIGHT  END  POINT  OF  THE  INTERVAL 
THE  NEXT  POINT  IN  THE  PROCESS.  IF  THERE  ARE 
NC  MCRE  POINTS  IN  THE  INTERVAL  (XLAST.ER). 
XNEXT  IS  ASSIGNED  THE  VALUE  ER  ♦ 1 .0  AND 
IER  IS  SET  AT  5. 

ERROR  FLAG.  IER  HAS  THE  FOLLOWING  MEANINGS: 

1.. .  LB  IS  NON-POSITIVE 

2. . .XLAST  IS  GREATER  THAN  ER 

3.. .GLB  IS  NEGATIVE 

5.. .  XNEXT  IS  GREATER  THAN  ER 


COMMENTS 

THE  INTENSITY  FUNCTION.  FCN,  IS  USER  SUPPLIED. 
EXAMPLE:  FUNCTION  FCN  (X  ) 

FCN  * 1.0  ♦ EXP(-X) 

RETURN 

.END 


PROGRAMM ER: 


LCDR  JOHN  SCOTT  REDO. 
AUG  1978 


USN 


SUBROUTINE  NHFNXT  ( I S . U e ,GLB. XLAST  , ER .XNEXT , I ER ) 

EXTERNAL  FCN 

CALL  CVFLCW 

CATA  EXPO/O. 0/, U/0. C/ 

RMIN  = GLB/UB 
IER  « 0 

IF  (UB.LE.G.O)  IER  - 1 
IF  (XLAST  .GE.ER)  I ER  = 2 
Ic  (GLB.LT.O.C)  I ER  *3 
IF  (IER.EC.l.OR.  IER.EQ.2)  RETURN 
IF  (IER. EC. 3)  GLB  » 0.0 


GENERATE  EMI.J  ); 
EXCEEDS  ER 


CHECK  TO  SEE  IF  ADDITION  OF  E*(I,J> 


EXMEAN  * 1.0/UB 
XNEW  * XLAST 


GENERATE  CNE  EXPONENTIAL  AND  SCALE 

i CALL  EXPON  (IS,  EXPO.l) 

EX  PG  * EX  PO*  EXMEAN 
XNEW  * XNEW+EXPO 
IF  (XNEW.GT.ER)  GO  TO  3 

GENERATE  UNIFORMIO,  1 ) THINNING  VARIATE 

CALL  RANDOM  (IS.U.l) 

TEST  LOWER  BOUND  FOR  THINNING 
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IF  (U.LE.RMIN)  GO  TO  2 

TEST  P0R  THINNING 

RATIO  = FCMXNEWJ/UB 
IF  (U.LE. RATIO)  GO  TO  2 
GO  TO  1 

ARRIVAL  HERE  INCICATES  SUCCESSFUL  THINNING 

2 XNEXT  * XNEto 
RETURN 


ARRIVAL  HERE  INDICATES  NO  MORE  POINTS  IN  INTERVAL 


3 IER  - 5 
XNEXT  * 
RETURN 
END 
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SUBROUTINE  DEGTWO 
SUBROUTINE  DEGTWO 


SHIS  PA8E  IS  BEST  QUALITY  PRACTICABL1 
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PURFCSE 

SI MULATE9  A NON-HOMOGENEOU S POISSON  PROCESS  WITH 
QUADRATIC  EXPONENTIAL  INTENSITY  FUNCTION  OVER 
A GIVEN  INTERVAL  USING  THE  POISSCN-DECOMPOS ITION 
AND  GAP  STATISTIC  ALGORITHM. 

USAGE 

CALL  DEGTWO (IS, A, Al, A2,EL,ER,II ,N, IER) 


DESCRIPTION  OF  PARAPETE  FS 
IS  - RANDOM  NUMBER  SEED. 

OR  LESS  DIGITS. 

A - CONSTANT  IN  INTENSITY 


INTEGER  WITH  NINE 


FUNCTION 


Al  - 1ST  DEGREE  COEFF  IN  INTENSITY  FUNCTION. 

A2  - 2ND  DEGREE  COEFF  IN  INTENSITY  FUNCTION. 

EL  - LEFT  END  POINT  OF  INTERVAL. 

ER  - RIGHT  END  POINT  OF  INTERVAL 
II  - 0 FOR  TIMES  OF  EVENTS. 

1 FOR  TIMES  BETWEEN  EVENTS. 

N - A VECTOR  OF  LENGTH  5.  N(l)  THROUGH  N(4) 
CONTAIN  NUMBERS  OF  EVENTS  FRCP  VARIOUS 
COMPONENTS  OF  THE  DECOMPOSED  INTENSITY 
FUNCTION.  N (5 ) CONTAINS  THE  TOTAL  NUMBER 
OF  EVENTS  IN  THE  NON-HOMOGENEOUS  POISSCN 
PROCESS  . 

COMMENTS 

CALLING  PRCGRAM  MUST  HAVE  A COMMON  REGION,  HOLD, 
OF  DIMENSION  (50C0),AND  AN  INTEGER  ARRAY  OF 
DIMENSION  (5). 

EXAMPLE:  DIMENSION  T(5000I,N(5) 

COMMON/HOLD/T 

CALLING  PROGRAM  MUST  CONTAIN  THE  FOLLWOING 
ASSIGNMENT  STATEMENT: 

M=N(5) 

CALLING  PROGRAM  MUST  USE  THE  FOLLOWING  JCL  CARDS 

//  EX  EC  FORT  CLG,  iMSL  = DP 
//FORT.SYSIN  DO  * 

TIMES  TO  EVENTS  OR  TIMES  BETWEEN  EVENTS  WILL  BE 
STORED  IN  CELLS  Till  THROUGH  T ( M)  . 

SUBROUTINE  DEGTWO  ( I S , A ,A1 , A2  ,EL ,E R , 1 1 , N, I ER ) 

Cl  MENS  ION  TTM  ES  ( 5 OOO  ),  T(5000),  N<5),  «><5> 

COMMON  /MIKE/  TIMES/HOLD/T 
CALL  OVFLOW 

INITIALIZE  VARIABLES 

P(l)  * A 
P( 2)  = Al 
P(  3 I = A2 
F(4)  = 0. 

P<  5)  * 0. 


DO  1 1=1,5 
1 N( I > = 0 


IF  RATE  FUNCTION  IS  LESS  THAN  DEGREE  TWO, 
USE  NHPP2  ROUTINE  ONLY 
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IF  (A2.EQ.0.)  GO  TO  2 
GO  TO  4 

CALL  NHPP2  (1$,EL,ER,A,A1,II,N1,IER.) 
N ( 5 ) * N1 

IF  (N1.EQ.0)  RETURN 


00  l 1=1,  N1 

tim£s<iT  = T ( I ) 
CONTINUE 


RETURN 

DETERMINE  COEFFICIENTS  FOR  MODIFIED 
DEGREE  ONE  RATE  FUNCTION 

4 TEST  = -A1/(2.*A2) 

TINT  = ER-EL 

IF  (A1.GE.0..AN0.A2.GT.  0. ) GO  TO  5 
GO  TO  6 

5 8 = A-A2*  TI  NT**2 
B1  = A1+2.*A2*TINT 
GO  TO  10 

6 9 = A 

IF  ( ( A1  .LE.O.  .AND.A2.LT. 0.1  .OR.  (A1.GT.0.  .AND  .A2  .LT.  0.  . 
1 TINT  ) ) GO  TO  7 
GO  TO  8 

7 81  = A 1+A  2*T  I NT 
GO  TO  10 

8 IF  (A1.GT.0. .AND.A2.LT. 0.. AND. TEST. LT. TINT!  GC  TO  9 
81  = Al 

GO  TO  10 

9 81  = Al/2. 

MUST  THE  INTERVAL  BE  PARTITIONED? 

10  IF  (Al*A2.LT.O.. AND.TEST.LT  .TINT)  C-0  TO  11 
ERNEW  * ER 

GO  TO  12 

11  ERNEW  = TEST+EL 

GENERATE  DEGREE  ONE  NHPP  ON  INTERVAL 

12  BB  « B 
BB1  * 01 

CALL  NHPP  2 (IStELfERNEWyBByBBlfOyNlflER) 

N(1  ) = N1 

IF  (N(l).EQkO)  GO  TO  14 

DO  13  I»1,N1 
TIMES! I)  * T ( I ) 

13  CONTINUE 

COMPUTE  LENGTH  OF  INTERVAL  AND  DETERMINE  VALUE 
OF  CSTAR  FOR  USE  IN  REJECTION  ROUTINE 


14  Q 


ii 

E3 

E4 

E5 

IF 

IF 

IF 

IF 

IF 


* ERNEW-EL 
= A 

= A2*0**2 
= A1*Q 

= Al **2/ (4 .*A2 ) 
= A1**2/(I2  .*A2  ) 


(Al.GE  .O..AND.A2.GT.O.)  GO  TO  15 
f Al.*  “ * “ " ‘ 


LT  .0.. AND. A2.GI.0.. AND. TEST. GE.TINI)  GC  TO  J6 


.TINT  ) GO  TO  19 


(Al.LT.O..  AND.  A2.GT.0.  .AND. TEST. LT. TINT)  GC  TO 
(Al.LE  .O..AND.  A2.LT.0.  ) GO  TO  18 
(AI.GT.0..AND.A2 .LT .0. .AND .TEST  .GE . 

CSTAR  = EXP<E1-E4)-EXP(  El) 

GO  TO  20 
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15  CSTAR  » EXP(iEi)— EXP  (El-  E2) 

GO  TO  20 

16  CSTAR  * EXP1E1  + E2  + E3  )-EXP(  E1  + E3 ) 

GO  TO  20 

17  CSTAR  = EXP(  E1-E4  )-EXP  ( E1-E5) 

GO  TO  20 

18  CSTAR  * EXP  ( E 1)-EXP  ( E1  + E3+E2) 

GO  TC  20 

19  CSTAR  » E XP( E1+  E3+  E2 )— EXP(E1 ) 

CGMPUTE  INTEGRAL  CF  MCCIFIEO  DEGREE  TWO  RATE  FUNCTION 
OVER  INTERVAL 

20  CALL  HELP  ( A,  A1  ,A2  , EL.ERNEW  .PMTR) 

PMTR  * PMTR-(EXP(B)*(EXP(B1*ERNEW)-EXP(B1*EL)  ))/Bi 

IDENTIFY  AS  FIRST  SUBINTERVAL 
NOTE  = 1 

GENERATE  REALIZATION  ON  POISON  (PMTR)  VARIATE 

21  CALL  PVAR  ( IS. PMTR, M) 

IF  (NCTE.EQ.l  ) GO  TO  22 
GO  TO  25 

REJECTION  ROUTINE  USEO  ON  FIRST  SUBINTERVAL 

22  N ( 2 ) * M 
P(  4 ) * B 
P(  5 I * B1 

IF  (N(2  ) • EQi.O  ) GO  TO  24 

CALL  REJECT  ( IS.EL.CSTAR.P, Q, N<2) ) 


00  23  I»1,M 
T IMES (N (1  U-ZI 
23  CONTINUE 


T ( I ) 


HAS  THE  INTERVAL  BEEN  PARTITIONED? 

24  IF  (ERNEW.EQ.ER)  GO  TO  34 
GO  TO  27 

USE  REJECTION  ROUTINE  ON  SECOND  PART  CF  INTERVAL 

25  N (4 ) * M 
P(  4 ) = B 
P<  5 ) * B 1 

IF  NO  EVENTS  OCCURRED  BYPASS  REJECTION  ROUTINE 

IF  (N(  4 ) • EQ'.O  ) GO  TO  35 
0 = ER-ELNEW 

CALL  REJECT  ( IS. ELNEW, CSTAR, P ,0. N< A ) ) 

COPY  TIMES  OF  EVENTS  INTO  ‘TIMES'  ARRAY 
N4  » N(l)+N(2)+N(3) 


DO  26  I *1 , M 
T IMES (N4+I ) * T(I) 
26  CONTINUE 


GENERATION  OF  VARIATES  COMPLETE. 
GO  TO  ORDERING  ROUTINE 

GO  TO  35 


oono  noon  non  noon  nnnn  on  noon  non 


IBIS  PASS  IS  BEST  QUALITY  FBJLCTICABLI 

um  (xxpy  puwistoD  to  ddc - 


c 

c 

c 

c 

c 

c 


INTERVAL  PARTITION  WAS  REQUIRED.  MUST  NOW 
CONSIDER  SECONO  SUBINTERVAL 

DETERMINE  COEFFICIENTS  FOR  MODIFIED 
DEGREE  ONE  RATE  FUNCTION 


27  IF  (A1.GT.0..AND.A2.LT.0.) 
E * A-A2*TINT**2 

B1  a A1*2.*A2*TINT 
GO  TO  29 

28  B = A+(A1/2.)*TINT 
B1  * Al/2.+A2*TI NT 

29  ELNEW  a ERNEW 

GENERATE  DEGREE 
B 


GO  TO  28 


ONE  NHPP  ON  INTERVAL 


BB 

B81  a B 1 

CALL  NJ-iPP2  CIS, ELNEW, ER, BB,  BB1,0,N2,IER) 
N(  3 ) a N3 

IF  (NOI.EQ.O)  GO  TO  31 
N3  * N(1 1 +Ni2  ) 

TRANSFER  TIMES  BETWEEN  ARRAYS 


DO  30  I-1,N3 
TIME$(N3  + I)  = T ( I ) 

30  CONTINUE 

31  Q » TINT 

determine  VALUE  OF  CSTAR 
THE  REJECTION  ROUTINE 


FOR  USE  IN 


E2  « A2*Q**2 
E3  » A1*Q 

IF  (A1.GT.0..AND.A2.LT.0.I  GO  TO  32 
CSTAR  « EXP(E1-E4)-EXP(E1-E5-E3-E2) 

GO  TC  33 

32  CSTAR  » EXP< E1-E4)-EXP(E1+E3+E2» 

COMPUTE  INTEGRAL  CF  MOCIFIED  OEGREE  TWO  RATE 
FUNCTION  OVER  SECOND  INTERVAL 

3?  CALL  HELP  ( A,  A1  , A2  , ELNEW  ,ER , PMT R ) 

PMTR  » PMTR-( EXPtB )*( E X P< B1*ER )- EXP < 81* ELNEW ) )) /B1 

IDENTIFY  AS  SECOND  SUBINTERVAL 

NOTE  = 2 
GO  TO  21 

PARTITION  OF  INTERVAL  NOT  REQUIRED.  COMPUTE  TOTAL 
EVENTS  AND  SUPERPOSE  TWO  EVENT  STREAMS 

34  N(5J  a N (1 ) +N (2 ) 

IF  (N(2).EQ.O)  GO  TO  38 
LBGN  « N( 1 ) +1 
JBGN  « 1 

CALL  COLATE  ( LBGN ,N < 5 ) , 1 ) 

GO  TO  38 

PARTITION  WAS  REQUIRED.  DETERMINE 
AMOUNT  OF  SORTING  NEEDED 

35  N (5  I * N(l)+N(2)+N(3)+N(4) 

IF  (N (2 ) . EQ. 0. AND. N (4) . EQ.O  i GO  TO  38 
IF  ( N ( 4 ) . EQ  • 0 ) GO  TO  36 
IF  (N(2  ) . EQLO  ) GO  TO  37 
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r * 
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MUST  SUPERPOSE  FOUR  EVENT  STREAMS 

L9GN  * N(  D + 1 
LFIN  * N!1)*N!2) 

CALL  COLATE  < LBGN , LFIN, 1 ) 

LBGN  = LF IN+N (3 )+l 
J3GN  * LF I N*1 
CALL  COLATE  I 
GO  TO  38 


MLST  SUPERPOSE  FIRST  HALF 

> 

LBGN, N2t 1 ) 


LBGN  »N  ( 5 ) y JBGN ) 

OF  ARRAY  ONLY 


36  N2  * NUI  + NI 
LBGN  = Nil* 
CALL  COLATE 
GO  TO  38 


+ N121 
mi 

TE  (I 


MLST  SUPERPOSE  SECOND  HALF  OF  ARRAY  ONLY 

37  KK  = N (1 ) +N  (2  )+l 

LBGN  = N!1)*N(2)+N!3)U 
LFIN  * N(  5) 

CALL  COLATE  I LBGN , N ( 5 ) , KK ) 

GO  TO  38 

ARE  TIMES  OF  EVENTS  OP  TIMES  BETWEEN  EVENTS  REOUESTED? 

38  IF  (II.EQ-04  RETURN 

CALCULATE  TIMES  BETWEEN  EVENTS 

S = TIMES  <11 
TIMES!!)  * TIMES!!  )— EL 
N5  = N( 5) 

CO  29  1=2 ,N5 
SI  * TIMES!  I ) 

TIMES!  I)  * TIMESU  )-S 
S = SI 

39  CONTINUE 

RETURN 

END 

SUBROUTINE  NHPP 2 SIMULATES  A NON-HO  MCGEN  EOUS 
PCISSCN  PROCESS  WITH  A LOG-LINEAR  INTENSITY 
IRATE)  FUNCTION 

SUBROUTINE  NHFP2  (IS»EL»ER,  A,  Al,  1 1 ,N , IER ) 

DIMENSION  T ( 5GOO) 

COMMON  /HOLD/  T 

CALL  OVFLOW 

INITIALIZE  VARIABLES 

IER  = 0 
TINT  * ER-EL 
A = EXP  I A +A1  *EL  ) 

IS  THE  POISSON  PROCESS  HOMOGENEOUS? 

IF  (A1.EQ.0.)  GO  TO  3 

PAR  = ! A* ( EXP l TINT* Al )-l« ) ) /A1 

IF  iAl.GT.O.)  GO  TO  1 

IFLAG  = 3 

GO  TO  2 

1 A = A*  EXP  (T  INT*A1  ) 

A1  = -A1 
IFLAG  = 2 
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C 
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COMPUTE  PARAMETERS  OF  BOTH  POISSON  RANCOM  VARIABLES 

2 THETA  * -A/Al 
GO  TO  4 

COMPUTE  RATE  AND  SCALING  PARAMETERS  FOR  HOMOGENEOUS 
PCISSON  PROCESS 

3 PAR  » TINT* A 
IpLAG  * l 
AANVRS  * I ./A 

COMPUTE  NUMBER  OF  EXPONENTIAL  VARIATES  REQUIRED 

4 NMA X = PAR+6.*SGRT(PAR) 

IS  THIS  A HOMOGENEOUS  POISSON  PROCESS? 

IF  (IFLAG.EQ.l)  GO  TO  17 
GENERATE  REALIZATION  ON  POISSON  (THETA)  VARIATE 

5 CONTINUE 

CALL  PVAR  ( IS.THETA,M) 

IF  (M.EQ.O)  GC  TO  7 

CALCULATE  TIMES  OF  EVENTS 

CALL  SEXPON  (IS,T»NMAX) 

B = -A1 
V * 0. 

JMA  > = NMAX+1 


I 

-■ 

1 

i 


DO  6 1=1 » JMA  X 

HAVE  NUMBER  OF  EVENTS  EXCEEDED  THE  MAXIMUM  NUMBER 
THAT  THE  ARRAY  CAN  HOLD? 

IF  (I.GT.NMAX)  GO  TO  8 
V = V+T  ( I )/  ( ( M— I+l ) *B  ) 

IF  ( V.GT .TINT  ) GO  TO  9 
TCI)  = V 

IF  (l.EQ.M)  GO  TO  10 
6 CONTINUE 


NO  EVENTS  OCCURRED 

7 N = 0 
RE  TURN 

TOO  MANY  EVENTS  FOR  ARRAY.  INCREMENT  ERROR 
COOE  AND  TRY  AGAIN 

8 IER  = IER+1 
GO  TO  5 

THE  NUMBER  OF  EVENTS  OBSERVED  TO  OCCUR  IN  THIS 
NON-HOMOGEN EOUS  POISSON  PROCESS  IS  ' N ' 

S N = 1-1 

IF  (N.EQ.O)  RETURN 
GO  TO  11 

10  N * M 

11  CONTINUE 

IS  THE  RATE  FUNCTION  INCREASING  OR  DECREASING? 

IF  CIFLAG.EQ.3)  GO  TO  13 

C TIME  REVERSAL  TECHNIQUE  IS  NECESSARY 
C DETERMINE  WHETHER  N IS  EVEN  OR  ODD 
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IS  I G * MOD  ( N»  2 ) 
NLOOP  = N/2 
N1  = N+i 


MIS  PACE  IS  BEST  QUALITY  PRACTICABLE 
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sa.ii,/rl’NL0CP 


T(  II  * ER-TIN1-I  ) 
T(NI-I)  * ER-S 
12  CONTINUE 


IF  (ISIG.EQ.l)  T(  NLCOP+1)  =ER— T (NLOQP+1) 
ARE  TIMES  OF  EVENTS  REQUESTED? 

FETURN 


IF  (II.EQ.O) 

GO  TO  15 
13  IF  (II.NE.O)  GO  TO  15 
IF  (EL.EQ.0,1  RETURN 


00  14  I *1 »N 
T( I)  = EL  +T ( I J 
14  CONTINUE 


RETURN 

CALCULATE  TIMES  BETWEEN  EVENTS 
15  S * Til) 


00  16  I =2  » N 
SI  * T(II 
Til)  = T( I)-S 
S = SI 
16  CONTINUE 


RETURN 

THE  PCISSON  PROCESS  IS  HOMOGENEOUS 


17  I = 1 
U = 0. 

CALL  SEXPCN  <IS,T,NMAX) 

18  L * U+T(I ) 

IF  (U.GT.PAR)  GO  TO  20 

i * m 

IF  (I.GT.NMAX)  GO  TO  19 
GO  TG  18 


INCREMENT  ERROR  CODE 
15  IER  = IER+1 

TRY  AGAIN  WITH  NEW  STRING  OF  VARIATES 


GO  TO  IT 
20  N * 1-1 

IF  (N  .EQ  .0)  RETURN 
IF  ( 1 1 « EQ  .1  I GO  TO  22 


DO  21  I =1  r N 
5L  = EL+AINVRS*T( I ) 
Till  * EL 
21  CONTINUE 


RETURN 
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22  00  23 
T ( I ) 

23  CONTINUE 


1=1  tN 

T(  I ) +AINVRS 


RETURN 

END 

SUBROUTINE  PVAR  GENERATES  A 
VARIATE,  M,  USING  THE  GAMMA 


POISSON 

METHOD 


(THETA) 


SUBROUTINE  PVAR  (IS, THETA, M) 
DIMENSION  T(5  000) 

COMMON  /HOLD/  T 
K = 0 
C * 16.0 
D = .875 

IF  (THETA. LT.C)  GO  TO  2 
GO  TO  5 
U = 1. 

CTN  = EXP (-THETA) 

MMA  X = THETA + 6.*  SORT  (THETA) 

I = 1 

CALL  SR AN D <IS,T,MMAX) 

U * U*T( I ) 

IF  (U.LT.CTN)  GO  TO  8 
I = 1 + 1 
K = K + l 

IF  (I.GT.MMAX)  GO  TO  3 
GO  TO  4 

NP  = INT ( 0*THET  A ) 

AN  = FLOAT ( NP ) 

CALL  GAMA  (AN,IS,G,1) 

IF  (G.GT. THETA)  GO  TO  6 

K = K+NP 

THETA  = THETA-G 

GO  TG  1 

L = THETA /G 

NP  = NP-i 

CALL  SRAND  (!IS,T,NP) 


DO  7 1=1, NP 
Ic  (T(I).LT.U) 
7 CONTINUE 


K = K+i 


C 

C 

C 

C 

c 


c 

c 


THE  VALUE  M IS  ASSUMED  BY  THE  POISSCN  (THETA)  VARIATE 

8 M = K 
RETURN 
END 

SUBROUTINE  REJECT  GENERATES  AN  ORDERED  SAMPLE 
OF  GIVEN  SIZE  FROM  THE  SECOND  COMPONENT 
OF  THE  ORIGINAL  INTENSITY  FUNCTION 
USING  A RE  JECTI  CN-ACCEPTANCE  ALGORITHM 


SUBROUTINE  REJECT 
DIMENSION  V( 5C0) ♦ 
DIMENSION  T ( 5G00) 
COMMON  /HCLD/  T 
L20  = L* 1 0 
IF  (L20.GT  .500)  L20 
LI  * L+l 
K = 1 
J = 0 

CALL  SRAND  (IS,V,L20) 


(IS,  EL,CSTAR,PVEC,Q,L  ) 


500 


DO 

J * 


2 1=1,  L2  0 
' J + l 
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C 

cc 


c 

c 


c 

c 


c 

c 


c 

c 

c 


TIKI  - Q*V(J)+EL 
J ■ J+i 

IF  (V  ( J ) .LT .CALC ( PVEC»  T (K ) ) /CSTAR ) K=K+1 
IF  (K.EQ.Ll)  GO  TO  3 
IF  ( J .GE .120- 1 ) GO  TO  I 

2 CONTINUE 

IF  (K.LT.L)  GC  TO  1 

3 CALL  PXSORT  <T,1,L> 

RETURN 

ENO 

SUBROUTINE  COLATE  SUPERPOSES  TWO  ORDERED 
EVE.IT  STREAMS  OVER  THE  SAME  INTERVAL 

SUBROUTINE  COLATE  1 LBGN  ,LFI  N , JBGN ) 
OiMENSl  ON  TIMES  (5000  )»  TI5000) 

COMMON  /MIKE/  TIME  S/HOLD/T 
I » JBGN 
J * I 
K * LBGN 

1 IF  (TIMESU  >.LT.TIMES(K  ))  GO  TO  2 
T(  J)  = TIMESIK) 

J * J+l 
K * K+l 

IF  (K.GT.LFIN)  GO  TO  3 
GO  TO  1 

2 T(  J i * TIMES  (I) 

J * J+l 

I = 1+1 
IF  (I.Ei 
GO  TO  1 

3 II  = LBGN-1 


LBGN)  GO  TO  5 


C 

c 


CO  4 N*I»  II 
T(  J ) = TIMES(N) 
J * J+l 

4 CONTINUE 

RETURN 

5 CONTINUE 

00  6 N=K»  LFI N 
TIN)  = TIMES(N) 

6 CONTINUE 


RETURN 

ENO 

SUBROUTINE  HELP  EVALUATES  THE  INTEGRATED  INTENSITY 
FUNCTION  OVER  TFE  INTERVAL  (EL,ER) 

SUBROUTINE  HELP  < A . A1 , A2»EL  .E R,  XX ) 

DOUBLE  PRECISION  MMCAW  » EB»  AA 
Z * SORT(ABS( A2) ) 

Y = <A1*Z  )/ <2.*A2> 

AA  » Z*E  L+Y 
SB  = Z*ER+Y 
CC  = A-A1*A1/(4.*A2) 

CC  » EXP( CC ) / Z 
IF  (A2.LT.0.)  GO  TO  1 

01  * CEXP  ( AA**2  )*MMOAW ( AA  ) 

02  » DEXP(B8**2)*MM0AW(BB) 

XX  a CC*(  Q2-Q1) 

RETURN 

1 CC  « CC*. 8862269  „ _ 

XX  » CC*(  DERF  (BB)-DERF(AA) ) 

RETURN 

ENO 

FUNCTION  CALC  EVALUATES  THE  SECOND  COMPONENT  OF 
THE  DECOMPOSED  INTENSITY  FUNCTION 


96 


I 


c 

c 


FOR  ANY  INPUT  VALUE. 

x » pum-pt 

XX  a P(4)  +p  . _ . 

CALC  a EXP(X)-EXP  (XX) 

RETURN 

END 


PU>«-P<?)*AB$A+P(3)*ABSA**2 
‘ Pi  5 )*A8  SA 
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